Least-Squares Estimation of the Integer GPS Ambiguities 勉強メモ

Least-Squares Estimation of the Integer GPS Ambiguities

P.J.G. Teunissen

Delft Geodetic Computing Centre (LGR) Department of the Geodetic Engineering Delft University of Technology Thijsseweg 11, 2629 JA DELFT The Netherlands

In [89]:
%matplotlib inline

import numpy as np
import numpy.matlib
import matplotlib.pyplot as plt
from   matplotlib.pyplot import *
import pandas as pd

from mpl_toolkits.mplot3d import Axes3D
from IPython.core.display import HTML , display


vec = lambda x: x[:,np.newaxis]
normalize = lambda v : v/np.linalg.norm(v)
meshgrid = np.meshgrid
pi = np.pi
sqrt = np.sqrt

1．最小二乗法

1.1 定式化

最小二乗法を使うときは・・・

過決定問題を解きたいときである． 過決定問題とは，未知数よりも独立な観測方程式の数が多い場合に未知数の最尤推定解を求める問題である．

過決定問題の例として・・・

二次元平面上で3点が与えられた時，それらをそれっぽく通る直線を引くことを考える．

In [90]:
#datas
p0,p1,p2=[1.,1.],[2.,2.],[3.,4.]
p=np.array([p0,p1,p2])

#table display
Label=(' u_i ',' v_i ')
Data=pd.DataFrame(p,columns=Label)
display(Data)

#plot datas
figure()
plot(p[:,0],p[:,1],'o')
xlim([0,5])
ylim([0,5])
grid('on')
xlabel('u')
ylabel('v')
u_i v_i
0 1.0 1.0
1 2.0 2.0
2 3.0 4.0
Out[90]:
Text(0, 0.5, 'v')

最小二乗法で必要なものは・・・

以下の3つである．

  • 未知数ベクトル
  • 観測値ベクトル
  • 未知数から観測値を予測するモデル(観測モデル,観測方程式)

今回，vを観測値として考える． また，観測モデルは直線関係とし，その傾きをα，切片をβとする．

  • 未知数ベクトル : $\boldsymbol{x}=[\alpha,\beta]^\mathbb{T}$
  • 観測値ベクトル : $\boldsymbol{y} = [v1,v2,v3]^\mathbb{T}$
  • 観測モデル : $v_\mathrm{i}=\alpha u_\mathrm{i}+ \beta$

観測モデルを行列表記すると・・・

観測ベクトルと未知数の関係を以下のように記述できる．(観測モデルを行列$\boldsymbol{A}$で表す) $$ \begin{eqnarray} \boldsymbol{y} &=& \boldsymbol{Ax} \\ \iff \left( \begin{array}{c} v1 \\ v2 \\ v3 \end{array} \right) &=&\left( \begin{array}{c} u1&1 \\ u2&1 \\ u3&1 \end{array} \right) \cdot \left( \begin{array}{c} \alpha \\ \beta \end{array} \right) \end{eqnarray} $$

この式を観測方程式と呼ぶ．

※未知数ベクトルと観測値ベクトルの関係が$\boldsymbol{y}=\boldsymbol{Ax}$ではなく$\boldsymbol{y}=\boldsymbol{A(x)}$のようにしか書けないとき，観測モデルは非線形である．非線形最小二乗法を使う必要がある．

In [91]:
A=np.c_[vec(p[:,0]),np.ones([3,1])]
y=vec(p[:,1])
print('A=' , A)
print('y=' , y)

A= [[1. 1.]
 [2. 1.]
 [3. 1.]]
y= [[1.]
 [2.]
 [4.]]

行列Aが正方行列でないので・・・

逆行列を求めることができない．なので，$\boldsymbol{x} = \boldsymbol{A}^{-1}\boldsymbol{y}$として解を求めたいが，それはできない． これは，以下の連立方程式を満たすα，βが存在しないことに対応する(多分) $$ \begin{eqnarray} v_\mathrm{1}&=&\alpha u_\mathrm{1}+ \beta \\ v_\mathrm{2}&=&\alpha u_\mathrm{2}+ \beta \\ v_\mathrm{3}&=&\alpha u_\mathrm{3}+ \beta \\ \end{eqnarray} $$

でも，各観測に誤差が乗っていると考えれば・・・

以下の連立方程式を満たすα，βを決定することができる． $$ \begin{eqnarray} v_\mathrm{1}&=&\alpha u_\mathrm{1}+ \beta + \varepsilon_1\\ v_\mathrm{2}&=&\alpha u_\mathrm{2}+ \beta + \varepsilon_2\\ v_\mathrm{3}&=&\alpha u_\mathrm{3}+ \beta + \varepsilon_3\\ \end{eqnarray} $$

ここで，$\varepsilon_i$は各観測に乗った誤差である．

誤差のある場合のモデルを行列形式で表示すると・・・

以下のようになる． $$\boldsymbol{y} = \boldsymbol{Ax} + \boldsymbol{\varepsilon}$$ ただし，$\boldsymbol{\varepsilon}$は以下の要素を持つ誤差ベクトルである． $$\boldsymbol{\varepsilon} = \left( \begin{array}{c} \varepsilon_1 \\ \varepsilon_2 \\ \varepsilon_3 \end{array} \right)$$

最小化二乗法の根本にある考え方は・・・

「誤差ベクトルはなるべく小さいほうが良いから，誤差ベクトルが最小になるように未知数ベクトルを選択しよう」というものである． 誤差ベクトルの大きさはノルムによって与えられる

ノルム (英: norm, 独: Norm) は、平面あるいは空間における幾何学的ベクトルの "長さ" の概念の一般化であり、ベクトル空間に対して「距離」を与えるための数学の道具である。(Wikipedia - https://ja.wikipedia.org/wiki/%E3%83%8E%E3%83%AB%E3%83%A0)

通常は，誤差の大きさ=誤差の二乗和=L2ノルムと考えて以下の式を最小化する．(別にL2ノルムでなくてもいいと思う．ノルムの公理を満たしていれば) $$||\boldsymbol{\varepsilon}|| = \sqrt{\varepsilon_1^2+\varepsilon_2^2+ \cdots + \varepsilon_n^2} = \boldsymbol{\varepsilon}^\mathrm{T} \boldsymbol{\varepsilon}$$

更にL2ノルムの概念を拡張して・・・

本項では，以下の式で表されるマハラノビス距離を採用する．マハラノビス距離を用いる最小二乗法を重み付き最小二乗法という． $$||\boldsymbol{\varepsilon}||_{\boldsymbol{Q}} = \sqrt{ \boldsymbol{\varepsilon}^\mathrm{T} Q^{-1} \boldsymbol{\varepsilon} }$$ ただし$$\boldsymbol{Q}=\mathrm{Cov}(\boldsymbol{y}) $$

マハラノビス距離は観測ベクトル$\boldsymbol{y}$の各要素同士の相関や分散を考慮した距離である．この距離を使うことで以下のようなシチュエーションに対応できる．

  • $v_1$に大きなノイズが乗っていることが予想されるから$v_1$をあまり信頼したくない
  • $v_2$のノイズは小さいと予想されるから$v_2$を信頼したい
  • $v_1$の値が大きいとき$v_3$の値も大きくなる傾向があるから，両者の関係を考慮したい
In [92]:
Qy=np.array([[5,0,3],[0,1,0],[3,0,1]])
print(Qy)

[[5 0 3]
 [0 1 0]
 [3 0 1]]

なお，$\boldsymbol{Q}=\mathrm{Cov}(\boldsymbol{x})=\boldsymbol{I}$(単位行列)のときマハラノビス距離はL2ノルム=ユークリッド距離に等しく，一般的な最小二乗法となる．

※「距離」は二点間の関係を表す量．ノルムはベクトルの大きさを表す量．誤差ベクトルの大きさを表すならマハラノビス距離ではなくマハラノビスノルムという言葉のほうが正しい気がする．しかしマハラノビスノルムという言葉は無いようだ．これ以降，誤差ベクトルの大きさを指してマハラノビス距離というときは，「誤差ベクトルの始点と終点間のマハラノビス距離」という意味合いで使う．

ここまでをまとめると重み付き最小二乗法は・・・

誤差ベクトルのマハラノビス距離を最小化する問題である．したがって，以下のように表される．

$$\mathrm{argmin} \ ||\boldsymbol{\varepsilon}||_\mathrm{Q} $$

ここで，$\boldsymbol{y} = \boldsymbol{Ax} + \boldsymbol{\varepsilon}$より，誤差ベクトルは以下で表される． $$ \boldsymbol{\varepsilon} = \boldsymbol{y} - \boldsymbol{Ax} $$

したがって，重み付き最小二乗法は以下のように定式化される． $$ \mathrm{argmin}_\mathrm{x} \ || \boldsymbol{y} - \boldsymbol{Ax} ||_\mathrm{Q}$$

マハラノビス距離は二乗根の演算が出てきて面倒くさいので二乗しておく $$ \mathrm{argmin}_\mathrm{x} \ || \boldsymbol{y} - \boldsymbol{Ax} ||^2_\mathrm{Q}$$

最小二乗法定式化 まとめ

Teunissen論文 2章 (1)式

線形重み付き最小二乗法は以下の最小化問題として定義できる． $$\mathrm{argmin}_\mathbf{x} ||\boldsymbol{y}-\boldsymbol{Ax}||^2_\mathrm{Q}$$ ただし $$\boldsymbol{x} \in \mathbb{R}^\mathrm{n}$$ $$\boldsymbol{y} \in \mathbb{R}^\mathrm{m}$$ $$m \geq n$$ $$\boldsymbol{A} \in \mathbb{R}^\mathrm{m \times n}$$ $$(\boldsymbol{A} \ : \ \mathbb{R}^\mathrm{n} \to \mathbb{R}^\mathrm{m})$$ $$\boldsymbol{Q}=\mathrm{Cov}(\boldsymbol{y}) \in \mathbb{R}^{\mathrm{m} \times \mathrm{m}}$$ $$||\cdot||^2_\mathrm{Q} = (\cdot)^\mathrm{T} \boldsymbol{Q}^{-1} (\cdot)$$

1.2 解法

表記の簡単のため，最小二乗問題の目的関数を以下のように記述する． $$J(\boldsymbol{x}) \equiv ||\boldsymbol{y}-\boldsymbol{Ax}||^2_\mathrm{Q} $$

また，目的関数を最小化するような未知数ベクトルを最尤推定解$\boldsymbol{\hat{x}}$と記述する． $$\boldsymbol{\hat{x}} \equiv \mathrm{argmin}_\mathbf{x} \ J(\boldsymbol{x})$$

In [93]:
J = lambda x,Q :  x.T@  np.linalg.solve(Q , x)

1.2.1 微分による解法

目的関数最小化のための一般的な考え方は・・・

「目的関数$J(\boldsymbol{x})$を未知数ベクトル$\boldsymbol{x}$で偏微分してゼロになるならば，$\boldsymbol{x}$は最尤推定解である」というものである． なぜなら，「ある未知数ベクトル$\boldsymbol{x}_0$が目的関数$J(\boldsymbol{x})$の最小値を与えるならば，$\boldsymbol{x_0}$ を微小変化させたときの目的関数$J(\boldsymbol{x})$の変化量がゼロになる」からだ． したがって，以下を満たす未知数ベクトルは，目的関数を最小化するための十分条件を満たす． $$\frac{\partial J(\boldsymbol{x})}{\partial \boldsymbol{x}} = 0$$

※微分してゼロって普通は目的関数の極大値や極小値や鞍点を与えるけれど，最小二乗法の場合はどうやら最小値を与えるらしい． (イメージ的には目的関数が誤差の二乗=下に凸の二次関数の形になる=極値は最小値ってことだと思う(未検証))． 実際に例題の目的関数値分布をプロットすると確かに極値=最小値っぽい．

※※今回の例題の目的関数値分布は最小値がたくさんあるように見えるけれど，拡大して見ると最小値が存在していることが分かる

In [94]:
#xi=[alpha,beta]^T : candidate of x
alp,bet = np.meshgrid(np.linspace(1,2,15),np.linspace(-0.6,0,15))#all of combination (matrix form)
xi=np.array([alp.ravel(),bet.ravel()])#all of combination (sets of 2x1 vector)
ei=np.matlib.repmat(y,1,alp.size) - A @ xi# resudual e : y-A*xi
Ji=np.diag(J(ei,np.eye(3)))#calc Costfunc for each xi : J(ei) @ Q=I
Jmesh=np.reshape(Ji,alp.shape)#Costfunc value (matrix form)

#plot J map
fig=figure()
ax = Axes3D(fig)
ax.plot_surface(alp,bet,Jmesh,cmap=plt.cm.coolwarm, rstride=1, cstride=1)#
ax.set_xlabel('alpha')
ax.set_ylabel('beta')
ax.set_zlabel('J(x)')
ax.view_init(30,110)

#zoom up (5% from the lowest)
idx=Jmesh<np.min(Jmesh)+(np.max(Jmesh)-np.min(Jmesh))*0.05
Az=np.linspace(125,90,2)
El=np.linspace(30,0,2)
for (az,el) in zip(Az,El):
    fig=figure()
    ax = Axes3D(fig)
    ax.scatter(alp[idx],bet[idx],Jmesh[idx])#
    ax.set_xlabel('alpha')
    ax.set_ylabel('beta')
    ax.set_zlabel('J(x)')
    ax.view_init(el,az)

目的関数$J(\boldsymbol{x})$がこのままでは微分しにくいので・・・

展開する．

$$ \begin{eqnarray} J(\boldsymbol{x})&=&||\boldsymbol{y}-\boldsymbol{Ax}||^2_\mathrm{Q} \\ &=& \left( \boldsymbol{y}-\boldsymbol{Ax} \right)^\mathrm{T} \boldsymbol{Q}^{-1} \left(\boldsymbol{y}-\boldsymbol{Ax} \right)\\ &=&\boldsymbol{y}^\mathrm{T} \boldsymbol{Q}^{-1} \boldsymbol{y} - \boldsymbol{y}^\mathrm{T} \boldsymbol{Q}^{-1} \boldsymbol{Ax} - (\boldsymbol{Ax})^\mathrm{T} \boldsymbol{Q}^{-1} \boldsymbol{y} + (\boldsymbol{Ax})^\mathrm{T} \boldsymbol{Q}^{-1} \boldsymbol{Ax} \\ &=&\boldsymbol{y}^\mathrm{T} \boldsymbol{Q}^{-1} \boldsymbol{y} - \boldsymbol{y}^\mathrm{T} \boldsymbol{Q}^{-1} \boldsymbol{Ax} - \boldsymbol{x}^\mathrm{T} \boldsymbol{A}^\mathrm{T} \boldsymbol{Q}^{-1} \boldsymbol{y} + \boldsymbol{x}^\mathrm{T} \boldsymbol{A}^\mathrm{T} \boldsymbol{Q}^{-1} \boldsymbol{Ax} \end{eqnarray} $$

目的関数$J(\boldsymbol{x})$を微分する

$$ \begin{eqnarray} &&\frac{\partial J(\boldsymbol{x})}{\partial \boldsymbol{x}} = 0 \\ &\iff& \frac{\partial}{\partial \boldsymbol{x}} \left( \boldsymbol{y}^\mathrm{T} \boldsymbol{Q}^{-1} \boldsymbol{y} - \boldsymbol{y}^\mathrm{T} \boldsymbol{Q}^{-1} \boldsymbol{Ax} - \boldsymbol{x}^\mathrm{T} \boldsymbol{A}^\mathrm{T} \boldsymbol{Q}^{-1} \boldsymbol{y} + \boldsymbol{x}^\mathrm{T} \boldsymbol{A}^\mathrm{T} \boldsymbol{Q}^{-1} \boldsymbol{Ax} \right) &=& 0 \\ &\iff& \left( 0 - \left( \boldsymbol{y}^\mathrm{T} \boldsymbol{Q}^{-1} \boldsymbol{A} \right)^\mathrm{T} - \boldsymbol{A}^\mathrm{T} \boldsymbol{Q}^{-1} \boldsymbol{y} + 2 \boldsymbol{A}^\mathrm{T} \boldsymbol{Q}^{-1} \boldsymbol{Ax} \right)　&=&0 \end{eqnarray} $$

ここで，ベクトルによる行列の微分に関する以下の公式を用いた．

$$ \begin{eqnarray} \frac{\partial (\boldsymbol{Ax}) } {\partial \boldsymbol{x}} &=& \boldsymbol{A}^\mathrm{T} \\ \frac{\partial (\boldsymbol{x}^\mathrm{T} \boldsymbol{A}) } {\partial \boldsymbol{x}} &=& \boldsymbol{A} \\ \frac{\partial (\boldsymbol{x}^\mathrm{T} \boldsymbol{Q} \boldsymbol{x} ) } {\partial \boldsymbol{x}} &=& \boldsymbol{Qx} \\ \end{eqnarray}$$

ただし，$\boldsymbol{A}$は一般の行列，$\boldsymbol{Q}$は対称行列を表す．

※微分前の左辺第4項の$\boldsymbol{A}^\mathrm{T} \boldsymbol{Q}^{-1} \boldsymbol{A}$は，対称行列である(証明別添)．

微分した式を整理する

ここで，共分散行列の逆行列$\boldsymbol{Q}^{-1}$は対称行列なので，$\boldsymbol{Q}^{-1}= \left( \boldsymbol{Q}^{-1} \right)^\mathrm{T} $であることを用いる．

$$ \begin{eqnarray} \left( 0 - \left( \boldsymbol{y}^\mathrm{T} \boldsymbol{Q}^{-1} \boldsymbol{A} \right)^\mathrm{T} - \boldsymbol{A}^\mathrm{T} \boldsymbol{Q}^{-1} \boldsymbol{y} + 2 \boldsymbol{A}^\mathrm{T} \boldsymbol{Q}^{-1} \boldsymbol{Ax} \right)　&=&0 \\ \iff 2\boldsymbol{A}^\mathrm{T} \boldsymbol{Q}^{-1} \boldsymbol{y} &=& 2 \boldsymbol{A}^\mathrm{T} \boldsymbol{Q}^{-1} \boldsymbol{Ax} \\ \iff \boldsymbol{x} &=& \left( \boldsymbol{A}^\mathrm{T} \boldsymbol{Q}^{-1} \boldsymbol{A} \right)^{-1} \boldsymbol{A}^\mathrm{T} \boldsymbol{Q}^{-1} \boldsymbol{y} \end{eqnarray} $$

したがって，最尤推定解$\boldsymbol{\hat{x}}$は以下の式で与えられる． $$\boldsymbol{\hat{x}} = \left( \boldsymbol{A}^\mathrm{T} \boldsymbol{Q}^{-1} \boldsymbol{A} \right)^{-1} \boldsymbol{A}^\mathrm{T} \boldsymbol{Q}^{-1} \boldsymbol{y}$$ また，観測に重みを付けない場合，$\boldsymbol{Q}=\boldsymbol{I}$より，最尤推定解$\boldsymbol{\hat{x}}$は以下の式で与えられる． $$\boldsymbol{\hat{x}} = \left( \boldsymbol{A}^\mathrm{T} \boldsymbol{A} \right)^{-1} \boldsymbol{A}^\mathrm{T}\boldsymbol{y}$$

In [95]:
xhat_W= np.linalg.solve((A.T @ np.linalg.solve(Qy,A)), A.T @ np.linalg.solve(Qy,y))
xhat=np.linalg.solve(A.T @ A, A.T @ y)
print('Normal LSE : xhat=\n',xhat)
print('-------------------------')
print('Weighted LSE : \n xhat=\n',xhat_W)
print('@ Weight Qy = \n',Qy)

Normal LSE : xhat=
 [[ 1.5       ]
 [-0.66666667]]
-------------------------
Weighted LSE : 
 xhat=
 [[ 1.625]
 [-1.125]]
@ Weight Qy = 
 [[5 0 3]
 [0 1 0]
 [3 0 1]]

以下に観測に重みを付けない場合の最小二乗直線を描画する．

In [96]:
t=np.arange(0,5,1)
figure()
plot(p[:,0],p[:,1],'o')
plot(t,xhat[0]*t+xhat[1],'r',label='LSE')
plot(t,xhat_W[0]*t+xhat_W[1],'b',label='WLSE')
xlim([0,5])
ylim([0,5])
grid('on')
xlabel('u')
ylabel('v')
legend()
Out[96]:
<matplotlib.legend.Legend at 0x1e98d261188>

LSEとWLSEを見比べる

WLSEに使った重み行列$\boldsymbol{Q}$を見ると，「v1のノイズが大きいから信頼したくない」というシチュエーションになっている． 実際に描画された最小二乗直線を見ると，確かに青いラインはv1を信頼していない(=最小二乗直線とv1の距離が大きくても気にしないことが)分かる

最小二乗法の微分解法 まとめ

最小二乗法の最尤推定解$\boldsymbol{\hat{x}}$は以下の式で与えられる． $$\boldsymbol{\hat{x}} = \left( \boldsymbol{A}^\mathrm{T} \boldsymbol{Q}^{-1} \boldsymbol{A} \right)^{-1} \boldsymbol{A}^\mathrm{T} \boldsymbol{Q}^{-1} \boldsymbol{y}$$

また，観測に重みを付けない場合，$\boldsymbol{Q}=\boldsymbol{I}$より，最尤推定解$\boldsymbol{\hat{x}}$は以下の式で与えられる． $$\boldsymbol{\hat{x}} = \left( \boldsymbol{A}^\mathrm{T} \boldsymbol{A} \right)^{-1} \boldsymbol{A}^\mathrm{T}\boldsymbol{y}$$

付録

$\boldsymbol{A}^\mathrm{T} \boldsymbol{Q}^{-1} \boldsymbol{A}$は，対称行列である．ただし，$\boldsymbol{Q}$は共分散行列

[証明]

共分散行列$\boldsymbol{Q}$は定義から実対称行列である． さらに，実対称行列の逆行列もまた実対称行列であるから，共分散行列の逆行列$\boldsymbol{Q}^{-1}$もまた実対称行列である．

ここで，実対称行列は下三角行列とその転置の積に分解できる(コレスキー分解)．したがって，ある上三角行列$\boldsymbol{U}$によって，共分散行列の逆行列は$\boldsymbol{Q}^{-1} = \boldsymbol{U}^\mathrm{T} \boldsymbol{U}$と表される． したがって左辺第4項は$\boldsymbol{A}^\mathrm{T} \boldsymbol{U}^\mathrm{T} \boldsymbol{U} \boldsymbol{A} = (\boldsymbol{UA})^\mathrm{T} \boldsymbol{UA}$である．

ここで，一般の行列$\boldsymbol{M}$について$\boldsymbol{M}^\mathrm{T} \boldsymbol{M}$は対称行列であるから，左辺第4項は対称行列である．

1.2.2 幾何学的解釈による解法

幾何学的とは具体的には・・・

m行n列の観測モデル$\boldsymbol{A}$を「n次元空間からm次元空間への線形写像」として見てみることである．

ここで，写像前後のそれぞれの空間に便宜的に以下のように名前を付ける．

  • 未知数空間 : 未知数ベクトルの存在するn次元空間
  • 観測空間 ： 観測ベクトルの存在するm次元空間

このとき，以下のことが成り立つ

  • 未知数空間の基底ベクトルは，$\boldsymbol{A}$の列ベクトルへと写像される
  • 未知数空間全体は，$\boldsymbol{A}$の列ベクトルを基底とする超平面に写像される

3点から直線を決定する問題で考える

$$ \begin{eqnarray} \boldsymbol{y} &=& \boldsymbol{Ax} \\ \iff \left( \begin{array}{c} v1 \\ v2 \\ v3 \end{array} \right) &=& \left( \begin{array}{c} u1&1 \\ u2&1 \\ u3&1 \end{array} \right) \cdot \left( \begin{array}{c} \alpha \\ \beta \end{array} \right) \end{eqnarray} $$

このとき，以下のような幾何学的解釈ができる．

  • 観測モデル$\boldsymbol{A}$は2次元空間から3次元空間への線形写像である．
  • 観測空間は3次元空間で，$v1,v2,v3$軸がある．
  • 未知数空間は2次元空間であり，$\alpha$軸，$\beta$軸がある．
  • 観測空間へと写像された未知数空間は3次元の観測空間内で平面(2次元)をなす．
  • 観測空間内における未知数平面は，$[u1,u2,u3]^\mathrm{T}$，$[1,1,1]^\mathrm{T}$を基底として持つ

この関係を以下に図示する．

In [97]:
xs=np.eye(2) # 2 vectors in Unknown Variable Space (UVS)
W=A @ xs # mapping to Observation Space (OS)
a,b=W[:,0],W[:,1] # 2 venctors on Unknown Variable Plane (UVP)
b_orth=b- ((a @ b)/(a @ a))*a# (a,b_orth) = 2 orthogonal vectors on UVP
e1,e2=vec(normalize(a)),vec(normalize(b_orth))# (e1,e2) = orthogonal basis of UVP
yhat=A @ xhat # projection obs y => UVP (to determaine UVP plot center)

#plane data for wireframe plot
param1,param2=meshgrid(np.linspace(-1,1,5),np.linspace(-1,1,5))
Wi=\
np.matlib.repmat(vec(param1.ravel()).T,3,1)*np.matlib.repmat(e1,1,param1.size)\
+np.matlib.repmat(vec(param2.ravel()).T,3,1)*np.matlib.repmat(e2,1,param2.size)\
+np.matlib.repmat(yhat,1,param1.size)#  Wi = param1*e1 + param2 *e2 + yhat
tomesh = lambda x : np.reshape(x,param1.shape)
Wx_mesh,Wy_mesh,Wz_mesh=tomesh(Wi[0,:]),tomesh(Wi[1,:]),tomesh(Wi[2,:])

for az in np.linspace(0,40,3):
    fig=figure()
    ax = Axes3D(fig)
    ax.plot(y[0],y[1],y[2],'om',label='observation')#
    ax.plot_wireframe(Wx_mesh,Wy_mesh,Wz_mesh,label='Unknown Variables Space')
    ax.quiver(yhat[0],yhat[1],yhat[2],W[0,0],W[1,0],W[2,0],\
              color='m',pivot='tail',label='Colum(A)=basis of UVS')
    ax.quiver(yhat[0],yhat[1],yhat[2],W[0,1],W[1,1],W[2,1],color='m',pivot='tail')
    ax.set_xlabel('v1')
    ax.set_ylabel('v2')
    ax.set_zlabel('v3')
    title('Observations Space')
    ax.view_init(5,az)
    legend(loc=3)

幾何学的に考えると誤差ベクトルは・・・

未知数平面上の任意の点と，観測値を結ぶベクトルである． なぜなら，そのベクトルの大きさは「観測モデルと観測値がどれくらい離れているか=観測モデル誤差の大きさ」に対応する尺度といえるからである．

したがって，誤差ベクトル$\boldsymbol{\varepsilon}$は以下のように定式化できる． $$\boldsymbol{\varepsilon} \equiv \boldsymbol{y} - \boldsymbol{Ax}$$

また，ベクトルの大きさをL2ノルムとすると誤差ベクトルの大きさは以下のように定式化される．

$$||\boldsymbol{\varepsilon}|| = \sqrt{ \boldsymbol{\varepsilon}^\mathrm{T} \boldsymbol{\varepsilon} }$$

最小二乗法を幾何学的に解釈すると・・・

誤差ベクトルの大きさ$||\boldsymbol{\varepsilon}||$を最小とするような未知数平面上の点$\boldsymbol{A\hat{x}}$を求める問題と捉えられる．

(a) 重みなしLSEの場合

誤差ベクトル最小となる未知数平面の点$\boldsymbol{\hat{y}}$は・・・

「未知数平面と誤差ベクトルが直交する」という条件を使って求める事ができる．

実際に「誤差ベクトルと未知数平面は直交するとき，誤差ベクトルのノルムが最小である」様子を図示する．

図示にあたって，直交する様子が分かりやすいように，未知数平面を真横から見ることを考える． 真横と言っても色々な方向があるので，今回は便宜上，以下の軸を設定する．

  • 横軸：LSEとWLSEの最尤推定解$\boldsymbol{\hat{x}},\boldsymbol{\hat{x}}_\mathrm{W}$に対応する平面上の点$\boldsymbol{\hat{y}},\boldsymbol{\hat{y}}_\mathrm{W}$を結んだベクトル方向
  • 縦軸：未知数平面からの高さ方向

この軸によって描画される平面をまず図示しておく．

In [98]:
yhat_W=A @ xhat_W
e1_=normalize(yhat_W-yhat) # horizontal axis base vector
e2_=normalize(y-yhat) # vertical axis base vector 

rng=np.linalg.norm(y-yhat_W) # plot range
param1,param2 = meshgrid(np.linspace(-rng,rng,25),np.linspace(-rng,rng,25))
valpoints=np.matlib.repmat(vec(param1.ravel()).T,3,1)*np.matlib.repmat(e1_,1,param1.size)\
+np.matlib.repmat(vec(param2.ravel()).T,3,1)*np.matlib.repmat(e2_,1,param2.size)\
+np.matlib.repmat(y,1,param1.size)# value points = param1*e1_ + param2 *e2_ + y
tomesh = lambda x : np.reshape(x,param1.shape)
valpx,valpy,valpz=\
tomesh(valpoints[0,:]),tomesh(valpoints[1,:]),tomesh(valpoints[2,:]),

AZ=np.linspace(13,110,3)
EL=np.linspace(16,62,3)
for (az,el) in zip(AZ,EL):
    fig=figure()
    ax = Axes3D(fig)
    ax.plot(y[0],y[1],y[2],'om',label='observation')#
    ax.plot(yhat[0],yhat[1],yhat[2],'oy',label='yhat = A*xhat (LSE)')#
    ax.plot(yhat_W[0],yhat_W[1],yhat_W[2],'oc',label='yhat = A*xhat (WLSE)')#
    ax.plot_wireframe(Wx_mesh,Wy_mesh,Wz_mesh)
    Y=np.c_[y,yhat]
    Y_W=np.c_[y,yhat_W]
    ax.plot(Y[0,:],Y[1,:],Y[2,:],'k')
    ax.plot(Y_W[0,:],Y_W[1,:],Y_W[2,:],'k')
    ax.plot_wireframe(valpx,valpy,valpz,rstride=50,cstride=50,label='Focus Plane')
    ax.set_xlabel('v1')
    ax.set_ylabel('v2')
    ax.set_zlabel('v3')
    legend(loc=0)
    ax.view_init(el,az)
    ax.dist=8

未知数平面を真横から見た状態で

「誤差ベクトルと未知数平面は直交するとき，誤差ベクトルのノルムが最小である」様子を以下に図示する．

In [99]:
# L2 norm from obs data to value points 
LSEnorm=np.diag( J(valpoints - np.matlib.repmat(y,1,param1.size) ,np.eye(3)))
LSEnorm_mesh=tomesh(LSEnorm)

# Viewpoint change (OS => Focus Plane)
P=np.c_[e1_,e2_,np.cross(e1_.T,e2_.T).T] # Rotation matrix 
transform = lambda v : np.linalg.solve(P ,(v-yhat)) # Coordinate transformation
Ty=transform(y)
Tyhat=transform(yhat)
Tyhat_W=transform(yhat_W)
Tvalpoints=transform(valpoints)
Tvalpx,Tvalpy=tomesh(Tvalpoints[0,:]),tomesh(Tvalpoints[1,:])


figure()
plt.contour(Tvalpx,Tvalpy,LSEnorm_mesh,17)
plot(Ty[0],Ty[1],'om',label='obs')
plot(Tyhat[0],Tyhat[1],'oy',label='yhat')
plot(Tyhat_W[0],Tyhat_W[1],'oc',label='yhat_W')
plot([Tyhat[0],Ty[0]],[Tyhat[1],Ty[1]],'k',label='normal to UVP')
axis('equal')
axis('square')
ylim(-0.01,1)
grid('on')
legend()
xlabel('UVPlane')
ylabel('hight from UVPlane')
title('L2 norm = const. contour')
Out[99]:
Text(0.5, 1.0, 'L2 norm = const. contour')

未知数平面と誤差ベクトルが直交するということは・・・

誤差ベクトルは未知数平面を張る基底ベクトルの全てと直交することから，内積を使って以下が成り立つ．

$$ \begin{eqnarray} ||\boldsymbol{y} - \boldsymbol{A\hat{x}}|| &\to& \mathrm{min}\\ \iff \left< \boldsymbol{a}_i , \boldsymbol{y} - \boldsymbol{A\hat{x}} \right> &=&\boldsymbol{0} \\ \iff \boldsymbol{a}_i^\mathrm{T} (\boldsymbol{y} - \boldsymbol{A\hat{x}}) &=&\boldsymbol{0}\\ \iff \boldsymbol{A}^\mathrm{T} (\boldsymbol{y} - \boldsymbol{A\hat{x}}) &=&\boldsymbol{0}\\ \iff \boldsymbol{A}^\mathrm{T} \boldsymbol{y} - \boldsymbol{A}^\mathrm{T} \boldsymbol{A\hat{x}} &=&\boldsymbol{0}\\ \iff \boldsymbol{\hat{x}} &=& (\boldsymbol{A}^\mathrm{T}\boldsymbol{A})^{-1} \boldsymbol{A}^\mathrm{T} \boldsymbol{y} \end{eqnarray} $$

ただし，$\boldsymbol{A} = [\boldsymbol{a}_1,\boldsymbol{a}_2, \cdots \boldsymbol{a}_i \cdots ]^\mathrm{T}$

このようにして最尤推定解を求めることができる．

(b) 重みつきLSEの場合

ノルムが通常のL2ノルムではなくマハラノビス距離なので・・・

未知数平面と誤差ベクトルの直交条件を使うことができない．

その様子を先程と同様に平面を真横から見た図で示す．

In [100]:
#mahalanobis norm from obs data to value points
WLSEnorm=np.diag( J(valpoints - np.matlib.repmat(y,1,param1.size) ,Qy))
WLSEnorm_mesh=tomesh(WLSEnorm)

figure()
plt.contour(Tvalpx,Tvalpy,WLSEnorm_mesh,29)
plot(Ty[0],Ty[1],'om',label='obs')
plot(Tyhat[0],Tyhat[1],'oy',label='yhat')
plot(Tyhat_W[0],Tyhat_W[1],'oc',label='yhat_W')
plot([Tyhat_W[0],Ty[0]],[Tyhat_W[1],Ty[1]],'k',label='NOT normal to UVP')
axis('equal')
axis('square')
ylim(-0.01,1)
grid('on')
legend()
xlabel('UVPlane')
ylabel('hight from UVPlane')
title('mahalanobis norm = const. contour')
Out[100]:
Text(0.5, 1.0, 'mahalanobis norm = const. contour')

幾何学的解釈によるLSE解法 まとめ

誤差ベクトルと未知数平面が直交する条件を用いるとLSEを解くことができる． しかし．WLSEでは直交条件が使えない．

1.2.3 幾何学的解釈の拡張とWLSE解法

前項ではWLSEでは誤差ベクトルと未知数平面との直交条件が使えないと述べた． ではWLSEの解はどのように導けばよいのか？

実は内積の概念を拡張すると・・・

WLSEでも直交条件が使えるようになる． 本項ではまず内積の拡張の準備として，線形空間，内積空間について述べる． 次に，WLSEでも直交条件が使えることを示すために，以下の定理を証明する．

「どんな内積空間においても， 任意の線形独立なベクトル$\boldsymbol{a,b}$に対して，その線形結合で表されるベクトル$\boldsymbol{a}-r\boldsymbol{b}$のノルムが最小となるような実数$r=r_0$を選ぶとき，ベクトル$\boldsymbol{a}-r_0\boldsymbol{b}$はベクトル$\boldsymbol{a}$と直交する．」

そして最後に，上記の定理を用いてWLSEでも直交条件が使えることを示す．

(a)ベクトル空間と内積空間

今まで扱ってきた「空間」というのは・・・

数学的には「標準内積空間」とか「ユークリッド内積空間」とかと呼ばれている． この空間は，我々が通常3次元空間と認識しているような空間で，ベクトルの大きさや2点間の距離がユークリッド距離を使って定義され，ベクトル間の角度がおなじみの内積を使って定義されている．

実は内積の定義によっては・・・

ベクトルの大きさがユークリッドノルムじゃないような内積空間も存在する． とにかく，内積が定義されていて，内積によって角度とか大きさが定義されている空間が一般的に内積空間と呼ばれている．

逆にいえば・・・

大きさとか角度とかが定義されていない空間というのも存在する．そういう空間はベクトル空間と呼ばれている． ではベクトル空間にはどんな概念が定義されているのか，というと，ベクトル同士の加法とベクトルのスカラー倍という演算だけ． きちんと数学的にベクトル空間を定義しようとするとアーベル群とか体とかいう言葉が出てきたりして自信がどんどんなくなっていくので，ざっくりとした数式にすると，少なくとも以下のような演算体系を持つ集合$V$がベクトル空間と呼ばれている． $$ \begin{eqnarray} \forall x,y \in V &,& \forall a\in \mathbb{R}\\ x+y &\in& V \\ ax &\in& V \end{eqnarray} $$

※メモ

ベクトル空間は上記の他にも分配法則とか結合法則とかいくつか成り立つべき条件があるっぽい． 加えて，集合Vに定義されている「+」という演算は逆元(1に対する-1)とか零元(0みたいに+をしても相手を変化させない数)とかが存在しなきゃいけないとか，交換法則が成り立っていなきゃいけないとか細かい条件がたくさんある． 更に更にスカラーの部分は実数じゃなくてもある一定の公理を満たす集合(おそらく体と呼ばれるもの)であればいいとかなんとか．．． なのでこのノートに書いてある議論は見る人が見たら不完全なはず．

身近なベクトル空間といえば・・・

複素数が挙げられる． なぜなら複素数同士の加法と複素数のスカラー倍もまた複素数だから． そうすると，「複素平面と二次元平面が図形的に対応する」みたいな感じに全く違う概念に対して共通の定理とか性質を適用することができるようになるから便利っぽい．

ベクトル空間のご利益って一体？ 線形写像とか，行列の固有値関係の定理とかはベクトル空間における議論だから，その上位の概念である内積空間においても必ず成立する．だからこれから内積の定義を拡張して遊んでも線形写像とか行列同士の積とかが今まで通り安心して使える．

(b)内積を用いたノルムと角度の定義

ベクトル空間に「内積」という演算を定義すれば内積空間になるけど・・・

さすがにどんな演算でも良いというわけではなく，以下の内積の公理を満たしている必要がある(公理の出どころは知らないけど)

$$ \forall x,y,z \in V , \forall a\in \mathbb{R}\\ \begin{eqnarray} \langle x,y \rangle &\in& \mathbb{R}\\ \langle x,x \rangle &\geq& 0 \\ \langle x,x \rangle =0 &\iff& x=0 \\ \langle x+y,z \rangle &=&\langle x,z \rangle +\langle y,z \rangle \\ \langle ax,y \rangle &=&a \langle x,y \rangle \\ \langle x,y \rangle &=&\langle y,x \rangle \\ \end{eqnarray} $$

上記の公理を満たす関数$\langle (\cdot),(\cdot) \rangle$を内積という． 概念的には，内積とは「何か二つのベクトルを一つの実数に対応付ける関数」なのだということが主な部分だと思う．

内積空間では，ノルムが楽に定義できる

三点$a,b,c$と，それぞれの点を結ぶベクトルを考える．また，ベクトルのノルムを$d(\cdot)$とする．このとき，ノルムの公理は以下のように表される．いずれも我々の直感に反さないようになっている．

  • 二点間の距離はマイナスではない．$ \iff d(\cdot)>=0$
  • 二点間の距離がゼロならば二点は一致する．$\iff \left[ d(\vec{ca})=0 \iff \vec{ca}=0 \right]$
  • a-b間の距離はb-a間の距離に等しい．$\iff d(\vec{ac})=d(\vec{ca})$
  • 回り道しないルートより回り道するルートのほうが距離は長い$ \iff d(\vec{ab})+d(\vec{bc}) \geq d(\vec{ac})$

ここで，$d(\cdot)=\sqrt{\langle \cdot, \cdot \rangle}$とするとなんとノルムの公理を満たす． したがって，$d(\cdot)=\sqrt{\langle \cdot, \cdot \rangle}$のように，内積を用いてノルムを定義することが可能である．表記の簡略化のため，$d(\cdot)$は通常$||\cdot||$と表される．

内積空間では，角度も楽に定義できる

数学の教科書で必ずといっていいほど載っている有名な「コーシー・シュワルツの不等式」によって，内積から角度が定義できる．

―Wikipedia 「コーシー・シュワルツの不等式」

x や y が実または複素内積空間 (X, 〈•, •〉) の元であるとき、 シュワルツの不等式は次のように述べられる： $$ |\langle x,y\rangle |^{2}\leq \langle x,x\rangle \cdot \langle y,y\rangle $$ 左辺は内積 〈x, y〉 の絶対値の平方である。

大切なのは，コーシー・シュワルツの不等式が任意の内積空間に対して成立するということ(証明略(というか分からない))． また，コーシー・シュワルツの不等式が後述のように角度を定義するから，それはつまり「ベクトル空間にどんな内積を定義しても角度が定義される」ということ． ベクトル空間に角度という概念を与える，なんていう根本的な定理だから有名なんだと思う．

角度の定義は，上述のノルムの定義を用いて，コーシー・シュワルツの不等式を以下のように変形することで得られる．(絶対値とノルムの記号が混在しているのでややこしい) $$ \begin{eqnarray} |\langle x,y\rangle |^{2} &\leq& \langle x,x\rangle \cdot \langle y,y\rangle \\ \iff |\langle x,y\rangle |^{2} &\leq& ||x||^2 ||y||^2 \\ \iff -||x|| \ ||y|| \leq &\langle x,y\rangle& \leq ||x|| \ ||y|| \\ \iff ||x|| \ ||y|| \ \cos \theta & =& \langle x,y\rangle \\ \iff \theta &=& \arccos \left( \frac{\langle x,y\rangle }{||x|| \ ||y|| \ } \right) \end{eqnarray} $$

この式見たことあるぞ！

In [101]:
def mdot(x,y,Q):
    return x.T@Q@y

def mnorm(x,Q):
    return sqrt(mdot(x,x,Q))

def mangle(x,y,Q):
    acos((mdot(x,y,Q))/(mnorm(x,Q)*mnorm(y,Q)))

ちなみに通常我々が使う内積は・・・

ユークリッド内積とか標準内積とか言われていて，以下の式で定義されている．

$$\langle \boldsymbol{x} , \boldsymbol{y} \rangle \equiv \boldsymbol{x} ^\mathrm{T} \boldsymbol{y} $$

なので，ノルムは以下のようなおなじみの式で定義することが可能である． $$ \begin{eqnarray} ||\boldsymbol{x}|| &\equiv& \sqrt{\boldsymbol{x}^\mathrm{T} \boldsymbol{x} } \\ &=& \sqrt{x_1^2 + x_2^2 + \cdots + x_n^2} \end{eqnarray} $$

(c) WLSEを考えるにあたって重要な内積空間の定理

重要な定理とはこれだ！

任意の内積空間において以下が成り立つ．

任意の線形独立なベクトル$\boldsymbol{a,b} \neq \boldsymbol{0} \in \mathbb{R}^n $に対して，その線形結合で表されるベクトル$\boldsymbol{a}-r\boldsymbol{b}$のノルムが最小となるような実数$r=r_0 \in \mathbb{R}$を選ぶとき，ベクトル$\boldsymbol{a}-r_0\boldsymbol{b}$はベクトル$\boldsymbol{a}$と直交する．つまり以下の式が成立する．

$$r_0=\mathrm{argmin}_\mathrm{r } ||\boldsymbol{b} -r \boldsymbol{a}|| \iff (\boldsymbol{b}-r_0 \boldsymbol{a}) \perp \boldsymbol{a} $$

この定理は内積の公理と角度の定義だけをつかって証明可能である．したがって，ユークリッド空間だけではなく，任意の内積空間に対して成り立つ．

[証明]

内積の公理を再掲する

$$ \forall x,y,z \in V , \forall a\in \mathbb{R}\\ \begin{eqnarray} \langle x,y \rangle &\in& \mathbb{R}\\ \langle x,x \rangle &\geq& 0 \\ \langle x,x \rangle =0 &\iff& x=0 \\ \langle x+y,z \rangle &=&\langle x,z \rangle +\langle y,z \rangle \\ \langle ax,y \rangle &=&a \langle x,y \rangle \\ \langle x,y \rangle &=&\langle y,x \rangle \\ \end{eqnarray} $$

まず，左辺を解く． $$r_0=\mathrm{argmin}_\mathrm{r } ||\boldsymbol{b} -r \boldsymbol{a}|| $$

目的関数を変形すると以下のようになる． $$ \begin{eqnarray} ||\boldsymbol{b} -r \boldsymbol{a}||&=& \sqrt{\langle \boldsymbol{b} -r \boldsymbol{a},\boldsymbol{b} -r \boldsymbol{a} \rangle} \\ &=&\sqrt{ \langle \boldsymbol{b},\boldsymbol{b} \rangle+ \langle \boldsymbol{b} , -r\boldsymbol{a} \rangle+ \langle -r\boldsymbol{a} ,\boldsymbol{b} \rangle + \langle -r\boldsymbol{a}, -r\boldsymbol{a} \rangle} \\ &=&\sqrt{ \langle \boldsymbol{b},\boldsymbol{b} \rangle -2r \langle \boldsymbol{b} , \boldsymbol{a} \rangle +r^2 \langle \boldsymbol{a} , \boldsymbol{a} \rangle } \end{eqnarray} $$

ここで，$\langle b,b \rangle \geq 0 $より，上記の式$f(r)$が最小になるということは以下の式が最小になることに等しい． $$ \begin{eqnarray} f(r)&=&-2r\langle \boldsymbol{b} , \boldsymbol{a} \rangle +r^2 \langle \boldsymbol{a} , \boldsymbol{a} \rangle \\ &=&\langle \boldsymbol{a} , \boldsymbol{a} \rangle \left( r^2 - 2r \frac{\langle \boldsymbol{b} , \boldsymbol{a} \rangle}{\langle \boldsymbol{a} , \boldsymbol{a} \rangle } \right)\\ &=&\langle \boldsymbol{a} , \boldsymbol{a} \rangle \left( r - \frac{\langle \boldsymbol{b} , \boldsymbol{a} \rangle}{\langle \boldsymbol{a} , \boldsymbol{a} \rangle } \right)^2 - \frac{\langle \boldsymbol{b} , \boldsymbol{a} \rangle ^2}{\langle \boldsymbol{a} , \boldsymbol{a} \rangle} \end{eqnarray} $$ したがって，目的関数が最小になるとき，すなわち，$r={\langle \boldsymbol{b} , \boldsymbol{a} \rangle}/{\langle \boldsymbol{a} , \boldsymbol{a} \rangle } = r_0$のとき以下が成立する． $$ \begin{eqnarray} f(r_0)&=&- \frac{\langle \boldsymbol{b} , \boldsymbol{a} \rangle ^2}{\langle \boldsymbol{a} , \boldsymbol{a} \rangle}\\ \iff -2r_0 \langle \boldsymbol{b} , \boldsymbol{a} \rangle +r_0^2 \langle \boldsymbol{a} , \boldsymbol{a} \rangle &=&- \frac{\langle \boldsymbol{b} , \boldsymbol{a} \rangle ^2}{\langle \boldsymbol{a} , \boldsymbol{a} \rangle} \\ \iff r_0^2 \langle \boldsymbol{a} , \boldsymbol{a} \rangle ^2 -2r_0 \langle \boldsymbol{b} , \boldsymbol{a} \rangle \langle \boldsymbol{a} , \boldsymbol{a} \rangle +\langle \boldsymbol{b} , \boldsymbol{a} \rangle ^2 &=&0\\ \iff \left( r_0\langle \boldsymbol{a} , \boldsymbol{a} \rangle - \langle \boldsymbol{b} , \boldsymbol{a} \rangle \right)^2 &=&0 \\ \iff \left( \langle \boldsymbol{a} , \boldsymbol{b} \rangle -\langle \boldsymbol{a} , r_0\boldsymbol{a} \rangle \right)^2 &=&0\\ \iff \langle \boldsymbol{a} , \boldsymbol{b} -r_0\boldsymbol{a} \rangle &=&0 \end{eqnarray} $$

ここで，2つのベクトル$\boldsymbol{a},\boldsymbol{b} -r_0 \boldsymbol{a}$のなす角を$\theta$とすると，角度の定義より以下が成り立つ． $$ \begin{eqnarray} \cos \theta &=& \left( \frac{\langle \boldsymbol{a},\boldsymbol{b} -r_0 \boldsymbol{a}\rangle }{||\boldsymbol{a}|| \ ||\boldsymbol{b} -r_0 \boldsymbol{a}|| \ } \right) \\ \iff \cos \theta &=& \left( \frac{0}{||\boldsymbol{a}|| \ ||\boldsymbol{b} -r_0 \boldsymbol{a}|| \ } \right)\\ \iff \theta &=& \frac{\pi}{2} \end{eqnarray} $$ したがって，以下の定理が任意の内積空間に対して成り立つ． $$r_0=\mathrm{argmin}_\mathrm{r } ||\boldsymbol{b} -r \boldsymbol{a}|| \iff (\boldsymbol{b}-r_0 \boldsymbol{a}) \perp \boldsymbol{a} $$

(d) WLSEを幾何学的に解く

WLSEってそもそも・・・

ユークリッドノルムではなく，マハラノビス距離(マハラノビスノルム)を最小にするような未知数を求める問題である．なので，ユークリッド内積空間での関数としてマハラノビス距離を考えるのではなく，マハラノビス距離がノルムの定義となるような内積空間を考えてみる．

マハラノビス距離がノルムの定義となるような内積空間を考えてみると・・・

内積の候補として以下の式が考えられる． $$\langle \boldsymbol{x},\boldsymbol{y} \rangle \equiv \boldsymbol{x}^\mathrm{T} \boldsymbol{Q}^{-1} \boldsymbol{y}$$ なぜなら，マハラノビス距離は以下の式で表されるからだ． $$||(\cdot)|| = \sqrt{ (\cdot)^\mathrm{T} \boldsymbol{Q}^{-1} (\cdot)}$$

そこで，上記の内積を仮にマハラノビス内積と呼び，マハラノビス内積によって定義される内積空間をマハラノビス内積空間と呼ぶ.

※どちらの言葉もグーグルに引っかからないので恐らく造語

満を持して・・・

誤差ベクトルのマハラノビス距離を最小にするWLSE問題を考える．

先程証明した定理より，いかなる内積空間においても，「誤差ベクトルのノルムが最小のとき，誤差ベクトルは未知数平面を張る基底ベクトルの全てと直交する」という命題が成り立つ．したがって以下の式が成り立つ．

$$ \begin{eqnarray} ||\boldsymbol{y} - \boldsymbol{A\hat{x}}||_\mathrm{Q} &\to& \mathrm{min}\\ \iff \boldsymbol{a}_i &\perp& (\boldsymbol{y} - \boldsymbol{A\hat{x}})\\ \iff \left< \boldsymbol{a}_i , \boldsymbol{y} - \boldsymbol{A\hat{x}} \right>_\mathrm{Q} &=&\boldsymbol{0} \\ \iff \boldsymbol{a}_i^\mathrm{T} \boldsymbol{Q}^{-1} (\boldsymbol{y} - \boldsymbol{A\hat{x}}) &=&\boldsymbol{0}\\ \iff \boldsymbol{A}^\mathrm{T} \boldsymbol{Q}^{-1} (\boldsymbol{y} - \boldsymbol{A\hat{x}}) &=&\boldsymbol{0}\\ \iff \boldsymbol{A}^\mathrm{T} \boldsymbol{Q}^{-1} \boldsymbol{y} - \boldsymbol{A}^\mathrm{T} \boldsymbol{Q}^{-1}\boldsymbol{A\hat{x}} &=&\boldsymbol{0}\\ \iff \boldsymbol{\hat{x}} &=& (\boldsymbol{A}^\mathrm{T}\boldsymbol{Q}^{-1}\boldsymbol{A})^{-1} \boldsymbol{A}^\mathrm{T} \boldsymbol{Q}^{-1}\boldsymbol{y} \end{eqnarray} $$

ただし，$\boldsymbol{A} = [\boldsymbol{a}_1,\boldsymbol{a}_2, \cdots \boldsymbol{a}_i \cdots ]^\mathrm{T}$

このようにして最尤推定解を求めることができる．

式がごちゃごちゃしてきた・・・

ので，整理する． 最尤推定解$\boldsymbol{ \hat{x}}$を観測モデル$\boldsymbol{A}$によって写像した点を$\boldsymbol{ \hat{y}}$とする．このとき，以下の式が成立する． $$\boldsymbol{ \hat{y}} = \boldsymbol{A \hat{x}} = \boldsymbol{A}(\boldsymbol{A}^\mathrm{T}\boldsymbol{Q}^{-1}\boldsymbol{A})^{-1} \boldsymbol{A}^\mathrm{T} \boldsymbol{Q}^{-1}\boldsymbol{y} $$

この式はよく見ると，観測ベクトル$\boldsymbol{y} $に適当な行列を掛け算することで$\boldsymbol{ \hat{y}} $へと写像しているように見える．ところで，$\boldsymbol{ \hat{y}} $は幾何学的に考えると観測ベクトルから観測モデル平面$\boldsymbol{Ax}$へと下ろした垂線との交点になっている． したがって，観測ベクトル$\boldsymbol{y} $を$\boldsymbol{ \hat{y}} $へと写像する行列は，平面$\boldsymbol{Ax}$への射影を行う行列であると言える． そこで，平面$\boldsymbol{Ax}$への射影行列を以下のように$\boldsymbol{P}_\mathrm{A}$と定義する． $$\boldsymbol{P}_\mathrm{A} \equiv \boldsymbol{A} \left( \boldsymbol{A}^\mathrm{T} \boldsymbol{Q}^{-1} \boldsymbol{A} \right)^{-1} \boldsymbol{A}^\mathrm{T} \boldsymbol{Q}^{-1} $$

このとき，$\boldsymbol{ \hat{y}} = \boldsymbol{P}_\mathrm{A} \boldsymbol{ y} $である． また，射影行列を用いると最尤推定解を以下のように表すことができる． $$\hat{\boldsymbol{x}} = \boldsymbol{A}^{-} \boldsymbol{P}_\mathrm{A} \ \boldsymbol{y}$$ ただし，$\boldsymbol{A}^{-}$は「$\boldsymbol{A}$の逆行列に近い概念を持つ適当な行列」で，以下が成り立つ． $$\boldsymbol{A}^{-}\boldsymbol{A}=1$$

また，射影行列を用いて，最小となる誤差ベクトルは以下のように表される． $$ \begin{eqnarray} \boldsymbol{\varepsilon} &=& \boldsymbol{y} - \boldsymbol{A \hat{x}} \\ &=&\boldsymbol{y} - \boldsymbol{P}_\mathrm{A} \boldsymbol{y} \\ &=&\left(\boldsymbol{I} - \boldsymbol{P}_\mathrm{A} \right)\boldsymbol{y} \\ &\equiv& \boldsymbol{P}_\mathrm{A}^\perp \boldsymbol{y} \end{eqnarray} $$

最終行の$\boldsymbol{P}_\mathrm{A}^\perp$は・・・

射影行列$\boldsymbol{P}_\mathrm{A}$に相補的な射影行列(正しい用法か怪しい)である．

どうやら，行列の性質を表す性質のひとつに「相補性」という言葉があるらしく，定義は足し算に関する逆行列のようなもの．適当な行列$\boldsymbol{A,B}$に対して$\boldsymbol{A}+\boldsymbol{B}=\boldsymbol{I}$が成立するとき，「$\boldsymbol{A}$は$\boldsymbol{B}$に対して相補的である」というらしい（A,B入れ替えても同じ）

そして，射影行列の定理の一つに，「射影行列に相補的な行列もまた射影行列である」というのがあるらしい．(参考 http://www.cis.twcu.ac.jp/~asakawa/waseda2002/math/linear-algebra.pdf @3.6_射影行列)

また，射影の性質に「直和分解可能である」というのがあるらしい．「(正確性に欠ける言い方だけれど，)あるベクトルを，別の二つの直交するベクトルの和として表現することが可能」という概念に近い性質だと思う．(Wikipedia_射影作用素)

※ここからは完全に予想になるけれど・・・

ある射影行列$\boldsymbol{P}$とそれに相補的な射影行列$\boldsymbol{P}^\perp$を使って適当なベクトル$\boldsymbol{x}$を各々射影すると，$\boldsymbol{Px} \perp \boldsymbol{P}^{\perp} \boldsymbol{x} $が成立するのだと思う．(なぜならばTeunissenがわざわざ$\perp$という記号を使っているから) 　

幾何学的解釈の拡張とWLSE解法 まとめ

Teunissen論文 2章(2)式

線形重み付き最小二乗法の解は以下で与えられる． $$\hat{\boldsymbol{x}} = \boldsymbol{A}^{-} \boldsymbol{P}_\mathrm{A} \ \boldsymbol{y}$$ ここで，$\boldsymbol{P}_\mathrm{A}$は観測ベクトルを観測モデル部分空間へ直交射影する行列であり，以下で表される． $$\boldsymbol{P}_\mathrm{A} = \boldsymbol{A} \left( \boldsymbol{A}^\mathrm{T} \boldsymbol{Q}^{-1} \boldsymbol{A} \right)^{-1} \boldsymbol{A}^\mathrm{T} \boldsymbol{Q}^{-1} $$ また，このとき以下が成り立つ $$\boldsymbol{\hat{y}}=\boldsymbol{P}_\mathrm{A} \ \boldsymbol{y}$$ $$\boldsymbol{\hat{\varepsilon}}=\boldsymbol{P}_\mathrm{A}^\perp \boldsymbol{y}$$ $$||\boldsymbol{\hat{\varepsilon}}||^2_{\mathrm{Q}} = ||\boldsymbol{P}_\mathrm{A}^\perp \boldsymbol{y}||^2_{\mathrm{Q}}$$ ただし $$\boldsymbol{P}_\mathrm{A}^\perp = \boldsymbol{I} - \boldsymbol{P}_\mathrm{A}$$

※ $\boldsymbol{A}^{-}$は「$\boldsymbol{A}$の逆行列に近い概念を持つ適当な行列」で，以下が成り立つ． $$\boldsymbol{A}^{-}\boldsymbol{A}=1$$ 実際の計算時には$\boldsymbol{A}^{-}$を直接求めることはせず，$\boldsymbol{A}^{-}\boldsymbol{A}=1$を利用して消去する．

2.整数最小二乗法

2.1 定式化

整数最小二乗法(ILS)とは

解が整数に限られるような状況下で最尤推定解を求める問題． ILSは以下のように定式化される．

$$\boldsymbol{\check{x}} = \mathrm{argmin}_{\mathrm{x}\in \mathbb{Z}} \ || \boldsymbol{y} - \boldsymbol{Ax} ||^2_\mathrm{Q}$$

ここで，$\boldsymbol{\check{x}}$は最尤推定整数解である．

ILSはNP困難問題

と，言われていて，効率的に解くのがなかなか難しいらしい．

解法の糸口として，ILS問題の最尤推定解$\boldsymbol{\check{x}}$は，LSE問題の最尤推定解$\boldsymbol{\hat{x}}$の近傍にあるということが使えそう．

なので，上記のILS式を$\boldsymbol{\hat{x}}$を含む式に変形することを考える．

そこで，準備として

任意の内積空間におけるピタゴラスの定理を証明する．

任意の内積空間において，$\boldsymbol{a}\perp\boldsymbol{b}$ならば，以下の式が成り立つ． $$||\boldsymbol{a}+\boldsymbol{b}||^2 = ||\boldsymbol{a}||^2 +||\boldsymbol{b}||^2$$

[証明]

$$ \begin{eqnarray} ||\boldsymbol{a}+\boldsymbol{b}||^2 &=& \langle\boldsymbol{a}+\boldsymbol{b},\boldsymbol{a}+\boldsymbol{b} \rangle\\ &=& \langle \boldsymbol{a},\boldsymbol{a} \rangle + 2\langle\boldsymbol{a},\boldsymbol{b} \rangle +\langle \boldsymbol{b},\boldsymbol{b} \rangle\\ &=& \langle \boldsymbol{a},\boldsymbol{a} \rangle +\langle \boldsymbol{b},\boldsymbol{b} \rangle \ \ \ \ \ (\because\ \boldsymbol{a}\perp\boldsymbol{b}) \\ &=&||\boldsymbol{a}||^2 +||\boldsymbol{b}||^2 \end{eqnarray} $$

ピタゴラスの定理を使うと

ILS問題の目的関数は以下のように変形できる．

$$\begin{eqnarray} || \boldsymbol{y} - \boldsymbol{Ax} ||^2_\mathrm{Q} &=& || \boldsymbol{P}_{\mathrm{A}} \boldsymbol{y} - \boldsymbol{Ax}||^2 + ||\boldsymbol{P}_\mathrm{A} \boldsymbol{y} - \boldsymbol{y}||^2\\ &=&|| \boldsymbol{A} \boldsymbol{\hat{x}} -\boldsymbol{A} \boldsymbol{x}||^2 + ||\boldsymbol{P}_\mathrm{A}^\perp \boldsymbol{y}||^2 \end{eqnarray}$$

ここで，右辺第２項はモデル化誤差ベクトルのノルムなので未知数ベクトル$\boldsymbol{x}$に依存しない．(観測モデルを設定した段階で確定する) したがって，ILS問題は以下のように定式化される

$$\begin{eqnarray} \boldsymbol{\check{x}} &=& \mathrm{argmin}_{\mathrm{x}\in \mathbb{Z}} \ || \boldsymbol{y} - \boldsymbol{Ax} ||^2_\mathrm{Q}\\ &=&\mathrm{argmin}_{\mathrm{x}\in \mathbb{Z}} \ || \boldsymbol{x} - \boldsymbol{\hat{x}} ||^2_\mathrm{Q} \end{eqnarray}$$

ここで，共分散行列$\boldsymbol{Q}$は任意の対称行列である．WLSE問題において採用した共分散行列は$\boldsymbol{Q}=\mathrm{Cov}(\boldsymbol{y})$だったことから，上記の変形後のILS問題では$\boldsymbol{Q}=\mathrm{Cov}(\boldsymbol{\hat{x}})$が適切であると考えられる．

なお，共分散行列$\boldsymbol{Q}=\mathrm{Cov}(\boldsymbol{\hat{x}})$は以下で与えられる．

$$\begin{eqnarray} \mathrm{Cov}(\boldsymbol{\hat{x}}) &=& \mathrm{Cov}\left( (\boldsymbol{A}^\mathrm{T}\boldsymbol{Q}^{-1}_\boldsymbol{y}\boldsymbol{A})^{-1} \boldsymbol{A}^\mathrm{T} \boldsymbol{Q}^{-1}_\boldsymbol{y} \boldsymbol{y}\right) \\ &=& (\boldsymbol{A}^\mathrm{T}\boldsymbol{Q}_\boldsymbol{y}^{-1}\boldsymbol{A})^{-1} \boldsymbol{A}^\mathrm{T} \boldsymbol{Q}^{-1}_\boldsymbol{y} \mathrm{Cov}(\boldsymbol{y}) \left( (\boldsymbol{A}^\mathrm{T}\boldsymbol{Q}^{-1}_\boldsymbol{y}\boldsymbol{A})^{-1} \boldsymbol{A}^\mathrm{T} \boldsymbol{Q}^{-1}_\boldsymbol{y} \right) ^\mathrm{T}\\ &=& (\boldsymbol{A}^\mathrm{T}\boldsymbol{Q}_\boldsymbol{y}^{-1}\boldsymbol{A})^{-1} \boldsymbol{A}^\mathrm{T} \boldsymbol{Q}^{-1}_\boldsymbol{y} \boldsymbol{Q}_\boldsymbol{y} \boldsymbol{Q}^{-1}_\boldsymbol{y} \boldsymbol{A} (\boldsymbol{A}^\mathrm{T}\boldsymbol{Q}^{-1}_\boldsymbol{y}\boldsymbol{A})^{-1} \\ &=& (\boldsymbol{A}^\mathrm{T} \boldsymbol{Q}^{-1}_\boldsymbol{y} \boldsymbol{A})^{-1} \end{eqnarray}$$

ここで，共分散に関する以下の性質を用いた． $$\mathrm{Cov}(\boldsymbol{Ax})=\boldsymbol{A}\mathrm{Cov}(\boldsymbol{x})\boldsymbol{A}^\mathrm{T}$$

整数最小二乗法 定式化 まとめ

第3章 式(3),(13)

ILS問題はLSE問題の最尤推定解$\boldsymbol{\hat{x}}$を用いて以下のように定式化される． $$\boldsymbol{\check{x}} = \mathrm{argmin}_{\mathrm{x}\in \mathbb{Z}} \ || \boldsymbol{x} - \boldsymbol{\hat{x}} ||^2_\mathrm{Q}$$ ここで，$\boldsymbol{\check{x}}$は整数最小二乗法の解である，

また，$\boldsymbol{Q}$は最尤推定解の共分散行列$\boldsymbol{Q}_\boldsymbol{\hat{x}}$であり，以下の式で与えられる． $$\boldsymbol{Q}_\boldsymbol{\hat{x}} = (\boldsymbol{A}^\mathrm{T} \boldsymbol{Q}^{-1}_\boldsymbol{y} \boldsymbol{A})^{-1}$$

※式(4)-(12)は省略． 未知数ベクトルに整数解と実数解が入り混じるとどうなるか(MLSE : Mixture LSE)，について述べている． 最終的にはMLSEはILS問題に帰着できるからILSさえ解ければ良いという結論のようだ．

2.2 解法

ILS問題の解き方は

基本的には解の候補抽出と評価の繰り返しである(しらみつぶし)． 言い換えれば，数式の変形の繰り返しによって代数的に解を求めることは困難である． ここがILSとLSEとの大きな違いである．

解候補抽出は

最初に探索範囲を決めておかないと，無限に抽出できてしまう． そこで，いま仮に目的関数が実数$\chi^2$以下になる領域を探索範囲として設定する． 数式的には以下のように表される．

$$|| \boldsymbol{x} - \boldsymbol{\hat{x}} ||^2_\mathrm{Q} \leq \chi^2$$

上記の式は，幾何学的には$\boldsymbol{\hat{x}}$を中心とする楕円の内側を表している． 原点中心のほうが今後の議論で扱いやすいため，以下のように変数変換を行う．

$$ \boldsymbol{q} \equiv \boldsymbol{x} - \boldsymbol{\hat{x}} $$

このとき，解候補の探索範囲は以下のように表される． $$ || \boldsymbol{q} ||^2_\mathrm{Q} \leq \chi^2 $$

当面の問題は，この楕円内にある整数解の候補をいかに効率よく列挙するか？である．

※ $ ||\boldsymbol{q}||=\mathrm{const.}$と言えば円の方程式である．なので，ノルムの定義を変えれば円も楕円も「とあるベクトルのノルム=const.」という単一の形で表現できる． これを利用すると，円について成り立つ様々な定理を「内積，角度，ノルム」を使って書き直して，内積をマハラノビス内積に置き換えることで，楕円について成り立つ定理が得られる． 例えば，楕円の接線を求めたいとき，まず円の接線方向ベクトル$\boldsymbol{a}$と「中心-接点」を結んだベクトル$\boldsymbol{b}$が直交することから，$\langle \boldsymbol{a} , \boldsymbol{b} \rangle =0$を得る．そして，内積をマハラノビス内積に置き換えると楕円の接線方程式が得られる．

In [102]:
def compose(Vinv,th): # get representation matrix
    R=np.array([[np.cos(th),-np.sin(th)],[np.sin(th),np.cos(th)]])
    return R@Vinv@R.T

def decompose(Qinv): # representation mat ->  Vinv,Rot(th)
    lam,R=np.linalg.eig(Qinv)
    Vinv=np.diag(lam)
    return (Vinv,R)

def abc2elip(a,b,chi2=1,resol=100,offset=np.array([[0],[0]])): #ax+by=chi^2
    th=np.linspace(-np.pi,np.pi,resol)
    x=np.c_[np.sqrt(chi2/a)*np.cos(th),np.sqrt(chi2/b)*np.sin(th)].T
    return (x+offset,np.diag([a,b]))

def abcth2elip(a,b,th,chi2=1,resol=100,offset=np.array([[0],[0]])):# th=radian
    x,Vinv=abc2elip(a,b,chi2,resol)
    R=np.array([[np.cos(th),-np.sin(th)],[np.sin(th),np.cos(th)]])
    return (R@x+ offset,Vinv,R,compose(Vinv,th))

def ivar2elip(Vinv,chi2=1,resol=100,offset=np.array([[0],[0]])): # (x^T)Qx=chi2  V=diag(var_x1,var_x2) //Quadratic form
    x,_=abc2elip(Vinv[0,0],Vinv[1,1],chi2,resol,offset)
    return x
    
def icov2elip(Qinv,chi2=1,resol=100,offset=np.array([[0],[0]])): # (x^T)Qx=chi2  Q=Cov(x)//Quadratic form
    Vinv,R=decompose(Qinv)
    x=ivar2elip(Vinv,chi2,resol)
    return (R@x+ offset , Vinv,R)


#elipsoid define
chi2=4 # radius
xhat=vec(np.array([1.1,2.1]))# LSE Solution
Qxinv=compose(np.array([[1/4,0],[0,1]]),pi/6) # Qxhat inv
Qx=np.linalg.inv(Qxinv)#Qxhat

#elipsoid draw
q,_,_=icov2elip(Qxinv,chi2=chi2)
x=q+np.matlib.repmat(xhat,1,q.shape[1])

figure()
plot(x[0,:],x[1,:])
plot(xhat[0,:],xhat[1,:],'ob',label='center=xhat with axis x')
grid('on')
legend()
xlabel('x1')
ylabel('x2')
figure()
plot(q[0,:],q[1,:])
plot(0,0,'ob',label='center = O with axis q')
grid('on')
xlabel('q1')
ylabel('q2')
legend()
Out[102]:
<matplotlib.legend.Legend at 0x1e98fb68488>

3. 補平面による候補点列挙

整数解候補の効率的列挙方法の一つに

補平面を利用するものがある． 補平面とは楕円を挟み込むような2つの平面である． もしそのような平面が見つかれば，その平面の間にある整数候補抽出は容易である． (補平面より解が-2.1から+1.4と分かれば解候補は-2,-1,0,1である)

※解候補抽出に必要な補平面数は2次元であれば2つの平面が2セット，n次元であれば2つの平面がnセットである．

3.1 接平面方程式

補平面を見つける方法の一つに

互いに平行な楕円の接平面を見つけるという方法がある． そのような接平面の方程式を考える．(参考: http://examist.jp/mathematics/space-equation/setuheimen/)

まず，円(球)への接平面方程式は，ベクトルを用いて以下の一般式で表される． $$(\boldsymbol{p}_0 -\boldsymbol{\hat{x}})^{\mathrm{T}}(\boldsymbol{x}-\boldsymbol{\hat{x}})=r^2$$ ここで，$\boldsymbol{p}_0$は円上の任意の位置を表すベクトル，$r$は定数である．

上記の式の左辺はユークリッド内積になっている．したがって，楕円の接平面方程式もアナロジーから同様に以下の式で表すことができると予想される．(証明していない) $$(\boldsymbol{p}_0 -\boldsymbol{\hat{x}})^{\mathrm{T}} \boldsymbol{Q}^{-1}(\boldsymbol{x}-\boldsymbol{\hat{x}})=\chi^2$$

ここで，$\boldsymbol{q}_0 = \boldsymbol{p}_0 -\boldsymbol{\hat{x}}$とおけば，$\boldsymbol{q}=\boldsymbol{x}-\boldsymbol{\hat{x}}$を用いて以下のように整理できる． $$\boldsymbol{q}_0^{\mathrm{T}}\boldsymbol{Q}^{-1}\boldsymbol{q}=\chi^2$$

まずはこの式が楕円の接平面を表すかどうか(力技で)確認する．

In [103]:
## support plane by vector on elipsoid
q0=vec(q[:,10]) #vector on elipsoid
print('q0 norm = chi =',mnorm(q0,Qxinv))

def classify(vectors,tfidx):
    '''
    vectors :NxM
    tfidx   :1xM bool
    return1 : NxM1 vectors (tfidx=true)
    return2 : NxM2 vectors (tfidx=false) (M1+M2=M)
    '''
    nidx=np.logical_not(tfidx)
    data_t=vectors[:,tfidx.ravel()]
    data_f=vectors[:,nidx.ravel()]
    return [data_t,data_f]

def mnorm_eval(Q,chi2):
    return lambda v1,v2 : mdot(v1,v2,Q) >= chi2

is_norm_gt_chi2=mnorm_eval(Qxinv,chi2)# function : norm greater than chi2 ? return bool

qq=(np.random.randn(2,400) )+np.matlib.repmat(q0,1,400)
qq_over,qq_under = classify(qq,is_norm_gt_chi2(q0,qq))

figure()
plot(q[0,:],q[1,:])
plot(qq_over[0,:],qq_over[1,:],'b.',label='greater than chi2')
plot(qq_under[0,:],qq_under[1,:],'.r',label='less than chi2')
plot([0,q0[0,:]],[0,q0[1,:]],'ok-',label='q0')
grid('on')
legend(loc=0)

q0 norm = chi = [[2.]]
Out[103]:
<matplotlib.legend.Legend at 0x1e98fc48888>

上図において，赤点と青点の境界が接線になっている

ように見える． そこで本ノートでは証明無しに以下の式が楕円の接線方程式であるとして話をすすめる． $$\boldsymbol{q}_0^{\mathrm{T}}\boldsymbol{Q}^{-1}\boldsymbol{q}=\chi^2$$

また，同じく図より，補平面から得られる解の候補は赤点で示す範囲だから，解の候補は以下の式で表される． $$\boldsymbol{q}_0^{\mathrm{T}}\boldsymbol{Q}^{-1}\boldsymbol{q} \leq \chi^2$$

より話を一般化するために(というかteunissenの論文上の式に寄せるために)

楕円上にないベクトルを用いた接線方程式を導出する． つまり，上記の接線方程式から，$\boldsymbol{q}_0$が楕円上の点であるという前提を取り除く．

いま，任意のベクトル$\boldsymbol{a}$を考える．ベクトル$\boldsymbol{a}$を用いて，楕円上のベクトル$\boldsymbol{q}_\mathrm{a}=f(\boldsymbol{a})$を作ることができれば，以下の楕円の接線方程式が得られる．

$$\boldsymbol{q}_{\mathrm{a}}^{\mathrm{T}}\boldsymbol{Q}^{-1}\boldsymbol{q}=\chi^2$$

楕円上のベクトルは

マハラノビスノルムが$\chi$である．したがって，以下の性質を満たすベクトルは楕円上のベクトルの一つである．

  • 方向単位ベクトル：$\frac{\boldsymbol{a}}{||\boldsymbol{a}||}$または$-\frac{\boldsymbol{a}}{||\boldsymbol{a}||}$
  • ノルム：$\chi$

上記の性質を満たすベクトルは以下で表される． $$\begin{eqnarray} \boldsymbol{q}_\mathrm{a} &=& \pm \chi \frac{\boldsymbol{a}}{||\boldsymbol{a}||}\\ &=&\pm \chi \frac{\boldsymbol{a}}{ \sqrt {\boldsymbol{a}^\mathrm{T} \boldsymbol{Q}^{-1} \boldsymbol{a} }} \end{eqnarray}$$

したがって，補平面方程式は

以下で表される．

$$\begin{eqnarray} \left( \pm \chi \frac{\boldsymbol{a}}{ \sqrt {\boldsymbol{a}^\mathrm{T} \boldsymbol{Q}^{-1} \boldsymbol{a} }} \right)_{\mathrm{a}}^{\mathrm{T}}\boldsymbol{Q}^{-1}\boldsymbol{q}&=&\chi^2 \\ \iff \boldsymbol{a}^{\mathrm{T}}\boldsymbol{Q}^{-1}\boldsymbol{q}&=& \pm \sqrt {\boldsymbol{a}^\mathrm{T} \boldsymbol{Q}^{-1} \boldsymbol{a} } \ \chi\\ \iff \boldsymbol{a}^{\mathrm{T}}\boldsymbol{Q}^{-1}(\boldsymbol{x}-\boldsymbol{\hat{x}})&=&\pm \sqrt {\boldsymbol{a}^\mathrm{T} \boldsymbol{Q}^{-1} \boldsymbol{a} } \ \chi \end{eqnarray}$$

また，補平面から得られる解の候補範囲は以下で表される． $$-\sqrt {\boldsymbol{a}^\mathrm{T} \boldsymbol{Q}^{-1} \boldsymbol{a} } \ \chi \leq \boldsymbol{a}^{\mathrm{T}}\boldsymbol{Q}^{-1}(\boldsymbol{x}-\boldsymbol{\hat{x}}) \leq \sqrt {\boldsymbol{a}^\mathrm{T} \boldsymbol{Q}^{-1} \boldsymbol{a} } \ \chi$$

In [104]:
######## support plane by any vector
a=np.random.rand(2,1)-0.5 #any vector
qa=sqrt(chi2)*a/mnorm(a,Qxinv) #vector on elipse

evalq=[]
for qai in [qa,-qa]:
    qq=(np.random.randn(2,800)) +np.matlib.repmat(qai,1,800)
    evalq.append(classify(qq,is_norm_gt_chi2(qai,qq)))

qq_over=np.c_[evalq[0][0],evalq[1][0]]
qq_under=np.c_[evalq[0][1],evalq[1][1]]

figure()
plot(q[0,:],q[1,:])
plot(qq_over[0,:],qq_over[1,:],'.b',label='greater than chi2')
plot(qq_under[0,:],qq_under[1,:],'.r',label='less than chi2')
plot([0,a[0,:]],[0,a[1,:]],'ok-',label='a')
plot([0,qa[0,:]],[0,qa[1,:]],'om-',label='qa')
plot([0,-qa[0,:]],[0,-qa[1,:]],'om-')
grid('on')
legend(loc=0)
Out[104]:
<matplotlib.legend.Legend at 0x1e98fbe1dc8>

3.2 法線ベクトルを用いた接平面方程式

補平面を用いた解候補抽出にあたって

接線が傾いていると候補抽出が面倒である． 逆に接線と軸が平行であれば，接線がy=-2.1，+1.4のような形になり，解候補が-2,-1,0,1であると簡単に分かる． そこで，接線と軸が平行(別の軸に対して直交)になるように任意のベクトル$\boldsymbol{a}$を選びたい．

「別の軸に対して直交」という条件は

接線の法線を考えれば簡単に適用できそうである．そこで，楕円上の法線ベクトルを考える．

一般的に，曲線$f(x_1,x_2)=0$上の点$f(X_1,X_2)$における法線ベクトルは$\left.(\frac{\partial f}{\partial x_1},\frac{\partial f}{\partial x_2})\right| _{x_1=X_1,x_2=X_2}$で与えられる．(http://mathtrain.jp/gradient) ベクトル形式で書くと，曲線$f(\boldsymbol{x})$上の点$f(\boldsymbol{x_0})$における法線は$\left. \frac{\partial f}{\partial \boldsymbol{x}} \right|_{\boldsymbol{x}=\boldsymbol{x_0}}$である．

したがって，楕円上の点$\boldsymbol{q}_\mathrm{a}$における法線ベクトル$\boldsymbol{c}$は以下で与えられる． $$\begin{eqnarray} \boldsymbol{c} &=& \left. \frac{\partial}{\partial \boldsymbol{x}} \left( \boldsymbol{x}^\mathrm{T} \boldsymbol{Q}^{-1} \boldsymbol{x} -\chi^2 \right) \right|_{\boldsymbol{x}=\boldsymbol{q}_\mathrm{a}} \\ &=&\left. 2\boldsymbol{Q}^{-1} \boldsymbol{x}\right|_{\boldsymbol{x}=\boldsymbol{q}_\mathrm{a}}\\ &=& 2\boldsymbol{Q}^{-1} \boldsymbol{q}_\mathrm{a}\\ &\propto& \boldsymbol{Q}^{-1} \boldsymbol{a} \end{eqnarray}$$

よって，任意のベクトル$\boldsymbol{c}$が法線ベクトルになるようなベクトル$\boldsymbol{a}$は以下で与えられる． $$\boldsymbol{a} =\boldsymbol{Qc}$$

楕円上の法線ベクトル式を接平面方程式に代入すると

任意のベクトル$\boldsymbol{c}$が法線ベクトルとなるような接平面方程式が以下のように得られる．

$$\begin{eqnarray} \boldsymbol{a}^{\mathrm{T}}\boldsymbol{Q}^{-1}(\boldsymbol{x}-\boldsymbol{\hat{x}})&=& \pm \sqrt {\boldsymbol{a}^\mathrm{T} \boldsymbol{Q}^{-1} \boldsymbol{a} } \ \chi \\ \iff \boldsymbol{c}^{\mathrm{T}} (\boldsymbol{x}-\boldsymbol{\hat{x}})&=& \pm \sqrt {\boldsymbol{c}^\mathrm{T} \boldsymbol{Q} \boldsymbol{c} } \ \chi \end{eqnarray}$$

接平面が軸に平行になるように取るには，軸方向単位ベクトルが法線になれば良い．そこで，i番目の軸方向単位ベクトルを以下の式で定義する． $$\boldsymbol{c}_\mathrm{i}=[0,0, \cdots 0,1,0, \cdots 0,0]^\mathrm{T}$$

実際に$\boldsymbol{c}_\mathrm{i}$を用いたときの補平面を以下に描画する． きちんと軸に沿った方向に接平面が生成されていることが分かる．

In [105]:
## support plane by normal vector
for column in np.eye(2):   
    c=vec(column)
    a=Qx@c
    qa=sqrt(chi2)*a/mnorm(a,Qxinv)

    evalq=[]
    for qai in [qa,-qa]:
        qq=(np.random.randn(2,800)) +np.matlib.repmat(qai,1,800)
        evalq.append(classify(qq,is_norm_gt_chi2(qai,qq)))
    
    qq_over=np.c_[evalq[0][0],evalq[1][0]]
    qq_under=np.c_[evalq[0][1],evalq[1][1]]
    
    
    figure()
    plot(q[0,:],q[1,:])
    plot(qq_over[0,:],qq_over[1,:],'.b',label='greater than chi2')
    plot(qq_under[0,:],qq_under[1,:],'.r',label='less than chi2')
    plot([0,a[0,:]],[0,a[1,:]],'ok-',label='a')
    plot([0,qa[0,:]],[0,qa[1,:]],'om-',label='qa')
    plot([0,-qa[0,:]],[0,-qa[1,:]],'om-')
    grid('on')
    legend(loc=0)

軸に平行な補平面方程式が得られたので

実際に解候補の抽出を行う．

単位ベクトルを接平面方程式に代入すれば,以下の式が得られる． $$ (x_\mathrm{i} - \hat{x}_\mathrm{i}) = \pm \sigma_\mathrm{\hat{x} \mathrm{i}} \chi$$

ただし，$\sigma^2_\mathrm{\hat{x} \mathrm{i}} = \boldsymbol{Q}_{\boldsymbol{\hat{x}}}[i,i]$である($\boldsymbol{Q}^{-1}$でないことに注意)

この式にしたがって，解候補抽出ができる．

In [106]:
lsub = lambda x: np.reshape(x,(1,x.size)) # x[:] matlab linear index

print('xhat=\n',xhat)
print('Qx=\n',Qx)
print('chi=\n',sqrt(chi2))
xmin,xmax=[],[] # tangent line
for i in [0,1]:
    xmin.append( -sqrt(chi2*Qx[i,i])+xhat[i])
    xmax.append( sqrt(chi2*Qx[i,i])+xhat[i])
    
xc_1=np.arange(np.ceil(xmin[0]),np.floor(xmax[0])+0.1)# x1 integer candidate
xc_2=np.arange(np.ceil(xmin[1]),np.floor(xmax[1])+0.1)# x2 integer candidate
print('x1_min=',xmin[0],',x1_max=',xmax[0])
print('x2_min=',xmin[1],',x2_max=',xmax[1])
print('xc_1=\n',xc_1)
print('xc_2=\n',xc_2)
xc_1m,xc_2m=meshgrid(xc_1,xc_2)
xc=np.r_[lsub(xc_1m),lsub(xc_2m)]  # candidate vector form

figure()
plot(x[0,:],x[1,:])
plot(xc[0,:],xc[1,:],'o',label='candidate')
grid()
legend()

xhat=
 [[1.1]
 [2.1]]
Qx=
 [[3.25       1.29903811]
 [1.29903811 1.75      ]]
chi=
 2.0
x1_min= [-2.50555128] ,x1_max= [4.70555128]
x2_min= [-0.54575131] ,x2_max= [4.74575131]
xc_1=
 [-2. -1.  0.  1.  2.  3.  4.]
xc_2=
 [-0.  1.  2.  3.  4.]
Out[106]:
<matplotlib.legend.Legend at 0x1e98d656348>

整数最小二乗法解法 補平面による候補点列挙 まとめ

第4章 式(15)(の下の文中の式)，(16)

ILS問題の目的関数が実数$\chi^2$以下になる領域を探索範囲として設定する． このときILS問題は，ひとまずは以下の式で与えられる楕円の内側にある整数ベクトル抽出問題に帰着される． $$|| \boldsymbol{x} - \boldsymbol{\hat{x}} ||^2_\mathrm{Q} \leq \chi^2$$ 整数解候補抽出の方法の一つに楕円内整数を挟むような互いに平行な補平面をを用いる方法がある．補平面方程式は以下で与えられる． $$\boldsymbol{c}^{\mathrm{T}} (\boldsymbol{x}-\boldsymbol{\hat{x}}) = \pm \sqrt {\boldsymbol{c}^\mathrm{T} \boldsymbol{Q} \boldsymbol{c} } \ \chi$$ ただし$\boldsymbol{c}$は任意のベクトルである．また，上記の補平面方程式の法線は$\boldsymbol{c}$である．補平面方程式に単位ベクトルを接平面方程式に代入すれば,以下の式で軸に沿った補平面が得られる． $$ (x_\mathrm{i} - \hat{x}_\mathrm{i}) = \pm \sigma_\mathrm{\hat{x}_i} \chi$$

4. 条件付き最小二乗法(CLS)を用いた候補点抽出

補平面による解候補抽出の特徴は

未知数ベクトルの要素を独立に求めることである．($x_1=[-2:4]$,$x_2=[0:4]$のように)

しかし，もしも解候補抽出の例において，$x_2=0$である，といった条件が仮に得られたとしたら，楕円方程式より，おおよそ$-2.3 \leq x_1\leq 1.4$という条件式が得られるはずだ． 本章では，未知数ベクトル同士の関係を考慮して，従属的に逐次的に解候補を求める方法について述べる．

In [107]:
##
D1=Qxinv[0,0]
L21=(Qxinv[1,0])/D1
D2=Qxinv[1,1]-(L21**2 *D1)

x2=0
q2=x2-xhat[1]
q1 =sqrt((D2/D1)*(chi2/D2 - q2**2)) - L21*q2
q1_=-sqrt((D2/D1)*(chi2/D2 - q2**2)) - L21*q2
x1,x1_=q1+xhat[0],q1_+xhat[0]
print(x1,x1_)

figure()
plot(x[0,:],x[1,:])
plot(xc[0,:],xc[1,:],'o',label='candidate')
plot([x1_,x1],[x2,x2],'o-',label='if x2 given')
grid()
legend()

[1.38040878] [-2.29810023]
Out[107]:
<matplotlib.legend.Legend at 0x1e98cff0c88>

4.1 二次元の場合

まずは簡略化のために

WLSE最尤推定解の共分散行列$\boldsymbol{Q}_{\boldsymbol{\hat{x}}}$の逆行列$\boldsymbol{Q}_{\boldsymbol{\hat{x}}}^{-1}$を以下の式で置き換える.

$$\boldsymbol{Q}_{\boldsymbol{\hat{x}}}^{-1} = \left( \begin{array}{c} a&b\\b&c\end{array} \right)$$

ただし，$a \gt 0, c \gt 0 ,ac-b^2 \gt 0$とする．

このとき，ILS問題の整数解候補は

以下の式を満たす． $$ \begin{eqnarray} \boldsymbol{q}^\mathrm{T}\boldsymbol{Q}_\boldsymbol{\hat{x}}^{-1}\boldsymbol{q} &\leq& \chi^2\\ \iff a q_1^2 + 2bq_1q_2 +cq_2^2 &\leq& \chi^2 \end{eqnarray}$$

ただし$q_\mathrm{i}=x_\mathrm{i} - \hat{x}_\mathrm{i}$である．

ここで，何らかの方法によって

$q_2$が与えられたと仮定すると，上記の式は$f(q_1)=\mathrm{const.}$の形に変形できる．具体的には以下のようになる．

$$\begin{eqnarray} a q_1^2 + 2bq_1q_2 +cq_2^2 &\leq& \chi^2 \\ \iff a \left(q_1^2 + \frac{2b}{a}q_1q_2 \right)+cq_2^2 &\leq& \chi^2 \\ \iff a \left(q_1+ \frac{b}{a}q_2\right)^2 -\frac{b^2}{a}q_2^2 +cq_2^2 &\leq& \chi^2 \\ \iff \left(q_1+ \frac{b}{a}q_2\right)^2 &\leq& \frac{1}{a}\left[ \chi^2 -\left(c -\frac{b^2}{a} \right)q_2^2 \right] \end{eqnarray}$$

どうやらこうやって最尤推定解を逐次絞り込んでいくやり方を条件付き最小二乗法(Conditional LSE)と呼ぶらしい(詳細不明)

上記の形式で与えられるCLSEについて

以下の定理が成り立つ．

$q_2=\mathrm{given}$の状態において与えられる$q_1$の範囲は3章で議論した解候補の範囲よりも少なくとも同等以下である．

[証明]

3章の補平面の議論より，解候補は以下の性質を満たす．

$$ (x_\mathrm{i} - \hat{x}_\mathrm{i})^2 \leq \sigma^2_\mathrm{\hat{x}_i} \chi^2$$

ただし，$\sigma^2_\mathrm{\hat{x} \mathrm{i}} = \boldsymbol{Q}_{\boldsymbol{\hat{x}} \mathrm{\_ i,i}}$である．ここで，$\boldsymbol{Q}_{\boldsymbol{\hat{x}}}$は以下の式で与えられる．

$$\begin{eqnarray} \boldsymbol{Q}_{\boldsymbol{\hat{x}}} &=& \left( \begin{array}{c} a&b\\b&c\end{array} \right)^{-1} \\ &=&\frac{1}{ac-b^2} \left( \begin{array}{c} c&-b\\-b&a\end{array} \right) \end{eqnarray}$$

したがって，補平面の考え方で得られる整数解候補は以下の範囲内にある整数である． $$\begin{eqnarray} q_1^2 &\leq& \frac{c}{ac-b^2}\chi^2 \\ q_2^2 &\leq& \frac{a}{ac-b^2}\chi^2 \end{eqnarray} $$

一方， $q_2=\mathrm{given}$の状態において$q_1$の整数解候補は以下範囲内にある整数である． $$\left(q_1+ \frac{b}{a}q_2\right)^2 \leq \frac{1}{a}\left[ \chi^2 -\left(c -\frac{b^2}{a} \right)q_2^2 \right]$$

ここで，「補平面の考え方で得られる$q_1$の探索範囲の大きさ」を$\alpha$で定義する．また，「$q_2=\mathrm{given}$の状態において得られる$q_1$の探索範囲の大きさ」を$\beta$と定義する．このとき，それぞれの探索範囲の大きさは以下で与えられる． $$\begin{eqnarray} \alpha &=& 2 \sqrt{ \frac{c}{ac-b^2}\chi^2 }\\ \beta &=& 2 \sqrt{\frac{1}{a}\left[ \chi^2 -\left(c -\frac{b^2}{a} \right)q_2^2 \right]}\\ \end{eqnarray} $$

上記の探索範囲の大きさについて，$\frac{\beta}{\alpha} \leq 1$ならば定理が成立する．ので確かめる．

$$\begin{eqnarray} \frac{\beta}{\alpha}&=&\frac{2 \sqrt{\frac{1}{a} \left[ \chi^2 -\left(c -\frac{b^2}{a} \right)q_2^2 \right]}}{2 \sqrt{ \frac{c}{ac-b^2}\chi^2 }}\\ &=&\frac{ \sqrt{\frac{1}{a} \left[ \chi^2 -\left(c -\frac{b^2}{a} \right)q_2^2 \right]}}{ \sqrt{ \frac{c}{ac-b^2}\chi^2 }}\\ &=&\frac{ \sqrt{ \color{red}{\frac{1}{a}} \chi^2 - \color{red}{\frac{1}{a}}\left(c -\frac{b^2}{a} \right)q_2^2 }}{ \sqrt{ \frac{c}{ac-b^2}\chi^2 }}\\ &=&\frac{ \sqrt{\frac{1}{a} \chi^2 - \frac{1}{a}\left(c -\frac{b^2}{a} \right)q_2^2 }}{ \sqrt{ \color{red}{\frac{1}{a-\frac{b^2}{c}}}\chi^2 }}\\ &\leq& \frac{ \sqrt{\frac{1}{a} \chi^2 - \frac{1}{a}\left(c -\frac{b^2}{a} \right)q_2^2 }}{ \sqrt{ \color{red}{\frac{1}{a}}\chi^2} }\\ &&\left(\because \frac{1}{a-\frac{b^2}{c}} \geq \frac{1}{a} \because \ \ \ b^2 \geq 0 , c \gt 0 \right)\\ &\leq& 1\\ &&\left( \because \frac{1}{a}\left(c -\frac{b^2}{a} \right)q_2^2 \geq 0 \ \ \ \because a \gt 0 , ac-b^2 \gt0 , q_2^2 \geq 0 \right) \end{eqnarray} $$

よって，$q_2=\mathrm{given}$の状態において与えられる$q_1$の範囲は3章で議論した解候補の範囲よりも少なくとも同等以下である．

上記の定理より

解候補の効率的な抽出にCLSEの考え方が有効であることが分かった．

4.2 多次元への拡張

2次元でのCLSEの考え方を多次元に拡張したい

2次元CLSEでは，未知数ベクトルの要素を一つずつ逐次を推定した． この考え方は連立方程式の後方代入という解法に似ている．(http://zakii.la.coocan.jp/signal/31_system_of_eqs.htm)

後方代入といえばコレスキー分解である

今回は分解したいものが対称行列なので，せっかくなので修正コレスキー分解(LDLt分解)する(結果ありきだけど) (https://en.wikipedia.org/wiki/Cholesky_decomposition#Computation)

任意の対称行列$\boldsymbol{A}$に対するLDLt分解:$\boldsymbol{A}=\boldsymbol{LDL}^\mathrm{T}$の一般式は以下で与えられる． $$\begin{eqnarray} \boldsymbol{D}_\mathrm{j} &=& \boldsymbol{A}_\mathrm{jj} - \displaystyle \sum_{k=1}^{j-1} \boldsymbol{L}^2_\mathrm{jk} \boldsymbol{D}_\mathrm{k} \\ \boldsymbol{L}_\mathrm{ij}&=&\frac{1}{\boldsymbol{D}_\mathrm{j}}\left( \boldsymbol{A}_\mathrm{ij}-\displaystyle \sum_{k=1}^{j-1} \boldsymbol{L}_\mathrm{ik}\boldsymbol{L}_\mathrm{jk}\boldsymbol{D}_\mathrm{k} \right)\ \ (\mathrm{for\ i \gt j}) \end{eqnarray}$$

In [108]:
# なんでnumpyにLDLt分解入ってないんだろう？
def LDLtdecomp(A):
    D,L=np.zeros(A.shape),np.eye(A.shape[0])
    for j in arange(0,A.shape[0]):
        D[j,j]=A[j,j]
        for k in arange(0,j):
            D[j,j]=D[j,j] - L[j,k]**2 *D[k,k]
        for i in arange(j+1,A.shape[0]):
            L[i,j]=A[i,j]/D[j,j]
            for k in arange(0,j):
                L[i,j]=L[i,j]-(L[i,k]*L[j,k]*D[k,k])/D[j,j]
    return (L,D)

二次元の場合，LDLt分解は以下のようになる

$$\begin{eqnarray} j=1 \\ \boldsymbol{D}_1 &=& \boldsymbol{A}_{11} - \displaystyle \sum_{k=1}^{0} \boldsymbol{L}^2_\mathrm{1k} \boldsymbol{D}_\mathrm{k}\\ &=&\boldsymbol{A}_{11}\\ j=1,i=2\\ \boldsymbol{L}_{21}&=&\frac{1}{\boldsymbol{D}_1}\left( \boldsymbol{A}_{21}-\displaystyle \sum_{k=1}^{0} \boldsymbol{L}_{2k}\boldsymbol{L}_{1k}\boldsymbol{D}_\mathrm{k} \right)\\ &=&\frac{\boldsymbol{A}_{21}}{\boldsymbol{D}_1}\\ j=2 \\ \boldsymbol{D}_2 &=& \boldsymbol{A}_{22} - \displaystyle \sum_{k=1}^{1} \boldsymbol{L}^2_\mathrm{2k} \boldsymbol{D}_\mathrm{k}\\ &=& \boldsymbol{A}_{22} - \boldsymbol{L}^2_{21} \boldsymbol{D}_1 \end{eqnarray}$$

楕円方程式にLDLt分解を適用してみる

今までの式を$\boldsymbol{A}=\boldsymbol{Q}_{\boldsymbol{\hat{x}}}^{-1}$と見れば，以下の式が得られる．

$$\begin{eqnarray} \boldsymbol{D}_1 &=& a \\ \boldsymbol{L}_{21} &=& \frac{b}{a}\\ \boldsymbol{D}_2 &=& c - \left( \frac{b}{a} \right) ^2 a\\ &=& c - \frac{b^2}{a} \\ \end{eqnarray}$$

よって楕円方程式は以下で与えられる． $$ \begin{eqnarray} \left( \begin{array}{c} q1\\q2\end{array} \right)^\mathrm{T} \boldsymbol{LDL}^\mathrm{T} \left( \begin{array}{c} q1\\q2\end{array} \right) &\leq& \chi^2\\ \iff \left( \begin{array}{c} q1\\q2\end{array} \right)^\mathrm{T} \left( \begin{array}{c} 1&0\\\frac{b}{a}&1 \end{array} \right) \left( \begin{array}{c} a&0\\0&c - \frac{b^2}{a} \end{array} \right) \left( \begin{array}{c} 1&0\\\frac{b}{a}&1 \end{array} \right)^\mathrm{T} \left( \begin{array}{c} q1\\q2 \end{array} \right) &\leq& \chi^2 \\a(q_1+\frac{b}{a}q_2)^2 + (c-\frac{b^2}{a})q_2^2&\leq& \chi^2 \end{eqnarray} $$

なんと都合よく先程見た形が出て来る．

なので，CLSEの多次元拡張はLDLt分解を使って表すことができそうである．

ここまでをまとめる

二次元楕円方程式は以下のように変形できる． $$\begin{eqnarray} \left( \begin{array}{c} q1\\q2\end{array} \right)^\mathrm{T} \left( \begin{array}{c} a&b\\b&c\end{array} \right) \left( \begin{array}{c} q1\\q2 \end{array} \right) &\leq& \chi^2 \\ \iff a(q_1+\frac{b}{a}q_2)^2 + (c-\frac{b^2}{a})q_2 &\leq& \chi^2 \end{eqnarray} $$

この変形はLDLt分解によって(恐らく高次元でも)導くことができる． $$ \begin{eqnarray} \left( \begin{array}{c} q1\\q2\end{array} \right)^\mathrm{T} \boldsymbol{LDL}^\mathrm{T} \left( \begin{array}{c} q1\\q2\end{array} \right) &\leq& \chi^2\\ \iff \left( \begin{array}{c} q1\\q2\end{array} \right)^\mathrm{T} \left( \begin{array}{c} 1&0\\\frac{b}{a}&1 \end{array} \right) \left( \begin{array}{c} a&0\\0&c - \frac{b^2}{a} \end{array} \right) \left( \begin{array}{c} 1&0\\\frac{b}{a}&1 \end{array} \right)^\mathrm{T} \left( \begin{array}{c} q1\\q2 \end{array} \right) &\leq& \chi^2\\ \iff a(q_1+\frac{b}{a}q_2)^2 + (c-\frac{b^2}{a})q_2^2&\leq& \chi^2 \end{eqnarray} $$

ここで，$q_2$が何らかの方法で与えられたと仮定すると，上記の式で得られる$q_1$の解候補範囲は，以下で与えられる． $$ \left(q_1+ \frac{b}{a}q_2\right)^2 \leq \frac{1}{a}\left[ \chi^2 -\left(c -\frac{b^2}{a} \right)q_2^2 \right]$$

この範囲は3章で議論した補平面で与えられる$q_1$の解候補の範囲よりも小さい．

これから考えること

高次元の場合もLDLt分解を用いれば，$q_{i+1},q_{i+2},\cdots,q_{n}$が与えられたとき，$q_{i}$の解候補の範囲を補平面の考え方よりも小さく取ることが可能ではないか？

アナロジーから多次元の場合を考える

楕円方程式はLDLt分解より，以下のように変形できる． $$\begin{eqnarray} \boldsymbol{q}^\mathrm{T}\boldsymbol{Q}_\boldsymbol{\hat{x}}^{-1}\boldsymbol{q} &\leq& \chi^2 \\ \iff \boldsymbol{q}^\mathrm{T}\boldsymbol{LDL}^\mathrm{T}\boldsymbol{q} &\leq& \chi^2 \\ \end{eqnarray} $$

ここで，$\boldsymbol{L}^\mathrm{T} \boldsymbol{q}=\boldsymbol{\tilde{q}}$とおく．このとき楕円方程式は以下で表される． $$\boldsymbol{\tilde{q}}^\mathrm{T}\boldsymbol{D}\boldsymbol{\tilde{q}} \leq \chi^2$$

更に，$\boldsymbol{D}=\boldsymbol{Q}_\mathrm{\boldsymbol{\tilde{q}}}^{-1}$とおけば，以下のように，あたかも未知数ベクトル$\boldsymbol{\tilde{q}}$の共分散行列が$\boldsymbol{Q}_\mathrm{\boldsymbol{\tilde{q}}}$であるかのように見える．

$$\boldsymbol{\tilde{q}}^\mathrm{T} \boldsymbol{Q}_\mathrm{\boldsymbol{\tilde{q}}}^{-1} \boldsymbol{\tilde{q}} \leq \chi^2$$

更に，$\boldsymbol{Q}_\mathrm{\boldsymbol{\tilde{q}}}^{-1}$は対角行列であるから，要素毎に書き下せば以下の式が得られる．

$$ \displaystyle \sum_{ i = 1 }^{ n } \frac{\boldsymbol{\tilde{q}}[i]^2}{\boldsymbol{Q}_\mathrm{\boldsymbol{\tilde{q}}}[i,i]} \leq \chi^2 $$

この式から，高次元の場合もLDLt分解を用いれば， $q_{i+1},q_{i+2},\cdots,q_{n}$が与えられたとき，$q_i$の解候補の範囲を補平面の考え方よりも小さく取ることが可能であると分かる．

例としてn=3のときを考える． $$ \frac{\boldsymbol{\tilde{q}}[1]^2}{\boldsymbol{Q}_\mathrm{\boldsymbol{\tilde{q}}}[1,1]} + \frac{\boldsymbol{\tilde{q}}[2]^2}{\boldsymbol{Q}_\mathrm{\boldsymbol{\tilde{q}}}[2,2]} + \frac{\boldsymbol{\tilde{q}}[3]^2}{\boldsymbol{Q}_\mathrm{\boldsymbol{\tilde{q}}}[3,3]} \leq \chi^2 $$

まず，$q_3$が仮に与えられた場合$\tilde{q}_3$は定数項になる．そこで，以下の式が与えられる． $$ \frac{\boldsymbol{\tilde{q}}[1]^2}{\boldsymbol{Q}_\mathrm{\boldsymbol{\tilde{q}}}[1,1]} + \frac{\boldsymbol{\tilde{q}}[2]^2}{\boldsymbol{Q}_\mathrm{\boldsymbol{\tilde{q}}}[2,2]} + \leq \chi^2 -\frac{\boldsymbol{\tilde{q}}[3]^2}{\boldsymbol{Q}_\mathrm{\boldsymbol{\tilde{q}}}[3,3]} $$

この式は結局楕円の式なので，(しかも半径が$\chi^2$よりも小さい)$$2次元楕円方程式の形に書くことができる． このとき，補平面の考え方を適用すれば，単純に補平面を適用するよりも小さい範囲が得られる． ここで，更に$q_2$を仮決定すれば左辺第二項が定数項になるので，以下略．

この記法を用いて楕円方程式を書き下せば，以下のようになる． $$\begin{eqnarray} \boldsymbol{\tilde{q}}^\mathrm{T} \boldsymbol{Q}_\mathrm{\boldsymbol{\tilde{q}}} \boldsymbol{\tilde{q}} &\leq& \chi^2 \\ \iff \displaystyle \sum_{i=1}^{n} \left( \frac{ x_\mathrm{i} - \hat{x}_{\mathrm{i \ | \ i+1,\cdots n}}}{\sigma_{\mathrm{i\ | \ i+1,\cdots,n}}} \right)^2 &\leq& \chi^2 \end{eqnarray}$$

ここで，変則的ではあるが， $$\boldsymbol{L}^\mathrm{T}=\left( \begin{array}{c} 1&(&&&&\boldsymbol{l}_1&) \\ 0&1&(&&&\boldsymbol{l}_2&) \\ 0&0&1&(&&\boldsymbol{l}_3&) \\ \vdots&&\ddots&\ddots&(&\boldsymbol{l}_i&)\\ \\ 0&0&\cdots&&&1 \\ \end{array} \right)$$

とおけば，$\boldsymbol{\tilde{q}} = \boldsymbol{L}^\mathrm{T} \boldsymbol{q}$は以下のように表される． $$ \boldsymbol{\tilde{q}} = \left( \begin{array}{c} q_1 + \boldsymbol{l}_1 [q_2,q_3,q_4\cdots,q_n]^\mathrm{T}\\ q_2 + \boldsymbol{l}_2 [q_3,q_4,\cdots,q_n]^\mathrm{T}\\ q_3 + \boldsymbol{l}_3 [q_4,\cdots,q_n]^\mathrm{T}\\ \vdots\\ q_{n-1} + \boldsymbol{l}_{n-1} [q_n]\\ q_n \end{array} \right) $$

なので，$q_{i+1},q_{i+2},\cdots,,q_{n}$が仮に与えられたときに，$\tilde{q}_i$に含まれる未知数は$q_{i}$のみになる．

CLSEの説明は終わったけれど更に式を整理したい

特に$\boldsymbol{\tilde{q}}$の要素がややこしいので以下の記法を導入する．

$$\begin{eqnarray} \tilde{q}_\mathrm{i} &=& q_i + \boldsymbol{l}_i [q_{i+1},\cdots,q_n]^\mathrm{T}\\ &\equiv & q_i + q_{\mathrm{i}\ | \ \mathrm{i+1,\cdots,n}} \end{eqnarray}$$

ここで，$q_i$は未知数，$q_{\mathrm{i}\ | \ \mathrm{i+1,\cdots,n}}$はgivenである．

更に，Teunissen論文に寄せるために，以下の記法を導入する．

1. $$\begin{eqnarray} \tilde{x}_\mathrm{i} &=& x_i - ( \hat{x} - q_{\mathrm{i}\ | \ \mathrm{i+1,\cdots,n}} )\\ &\equiv& x_i - \hat{x}_{\mathrm{i}\ | \ \mathrm{i+1,\cdots,n}} \end{eqnarray}$$

2. $$\boldsymbol{Q}_\mathrm{\boldsymbol{\tilde{q}}}[i,i] = \frac{1}{\boldsymbol{D}[i,i]} = \frac{1}{ \sigma^2_{\mathrm{i}\ | \ \mathrm{i+1,\cdots,n}}}$$

この記法を用いて楕円方程式を書き下せば，以下のようになる． $$\begin{eqnarray} \boldsymbol{\tilde{q}}^\mathrm{T} \boldsymbol{Q}_\mathrm{\boldsymbol{\tilde{q}}}^{-1} \boldsymbol{\tilde{q}} &\leq& \chi^2 \\ \iff \displaystyle \sum_{i=1}^{n} \left( \frac{ x_\mathrm{i} - \hat{x}_{\mathrm{i \ | \ i+1,\cdots n}}}{\sigma_{\mathrm{i\ | \ i+1,\cdots,n}}} \right)^2 &\leq& \chi^2 \end{eqnarray}$$

プログラムでCLSEを確認する

$x_2$の候補を補平面方程式から求め，CLSEを用いて$x_1$の候補を列挙する．

In [ ]:

LAMBDA法

LDLt変換 $$\boldsymbol{Q}^{-1}_{x}=\boldsymbol{LDL}^\mathrm{T}$$ より $$\boldsymbol{Q}_{x}=(\boldsymbol{L^\mathrm{T}})^{-1} \boldsymbol{D}^{-1}\boldsymbol{L}^{-1}$$

また，Z変換により $$\boldsymbol{Q}_{z}=\boldsymbol{Z}^\mathrm{T} \boldsymbol{Q}_{x} \boldsymbol{Z}$$

合わせると $$\boldsymbol{Q}_{z}=\boldsymbol{Z}^\mathrm{T}(\boldsymbol{L^\mathrm{T}})^{-1} \boldsymbol{D}^{-1}\boldsymbol{L}^{-1} \boldsymbol{Z}$$

なので， $$ \boldsymbol{Z}^\mathrm{T}(\boldsymbol{L^\mathrm{T}})^{-1} \approx \boldsymbol{I}\\ \boldsymbol{L}^{-1} \boldsymbol{Z}\approx \boldsymbol{I} $$ となるようなユニモジュラZを取れば変換後はほぼ無相関なデータ空間における議論になる．